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PRINCIPAL NOTATION 


SYMBOLS 


Unless specified otherwise, all variables are nondimensional. 


Symbol 


Definition 


A, B, C 
A', B', C' 

Cpy Cp 

E, F, G 


E t 



Ekj 

Fk, 

G 

A 

A 

A 

Ek, 

Fk, 

G 


AAA 


Ek, 

, Fy-p Gk[ 

Ek 2 

. Fy 2 , G V2 

F, 

G, H 

A 

F, 

G, H 


hr 

ij> k 
J 
k 
k 


k h k t 

Lr 

N eq 

N u N 2 , Nz 

P 

Pr r 

Pri , Pr t 


Speed of sound. 

Coefficient submatrices in block tridiagonal system of equations. 

Coefficient submatrices for boundary conditions. 

Specific heats at constant pressure and volume. 

Inviscid flux vectors in the Cartesian coordinate form of the governing equations. 

Inviscid flux vectors in the computational coordinate form of the governing 
equations. 

Total energy per unit volume. 

Viscous flux vectors in the Cartesian coordinate form of the governing equations. 

Viscous flux vectors in the computational coordinate form of the governing 
equations. 

Non-cross derivative viscous flux vectors in the computational coordinate form of 
the governing equations. 

Cross derivative viscous flux vectors in the computational coordinate form of the 
governing equations. 

Flux vectors in the Cartesian coordinate form c n the k-z turbulence model 
equations. 

Flux vectors in the computational coordinate form of the k-z turbulence model 
equations. 

Stagnation enthalpy per unit mass. 

Grid indices in the t], and £ directions. 

Jacobian matrix of the generalized grid transformation. 

Effective thermal conductivity coefficient. 

Turbulent kinetic energy. 

Laminar and turbulent thermal conductivity coefficient. 

Dimensional reference length. 

Number of governing equations being solved. 

Number of grid points in the £, rj, and £ directions. 

Static pressure. 

Reference Prandtl number. 

Laminar and turbulent Prandtl number. 

Heat fluxes in the Cartesian x, y , and z directions. 


Proteus 3-D Analysis Description 


Principal Notation 3 


Symbol 


Definition 


Q Vector of dependent variables in the Cartesian coordinate form of the governing 

equations. 

A Vector of dependent variables in the computational coordinate form of the gov- 

v eming equations. 

R Gas constant. 

Re r Reference Reynolds number. 

S Source term subvector in block tridiagonal system of equations. 


S' 

S, T 


S, T 

t 


Source term subvector for boundary conditions. 

Non-derivative terms in the Cartesian coordinate form of the k-z turbulence model 
equations. 

Non-derivative terms in the computational coordinate form of the k-z turbulence 
model equations. 

Physical time. 


T 

W, V, w 

W 


A 

w 


x,y, z 


Static temperature. 

Velocities in the Cartesian y , and z directions. 

Vector of dependent variables in the Cartesian coordinate form of the k-z turbu- 
lence model equations. 

Vector of dependent variables in the computational coordinate form of the k-z 
turbulence model equations. 

Cartesian coordinates. 


y 

6 

A, V 


eg>, 

£/ 

ef>, sf\ etc. 

Ou O3 


k 2 , k 4 





Mr 

v 




p 

a 


Ratio of specific heats, c^/c,. 

Difference operator. 

First-order forward and backward difference operators. 

Turbulent dissipation rate. 

Second- and fourth-order explicit artificial viscosity coefficients in constant coeffi- 
cient model. 

Implicit artificial viscosity coefficient. 

Second- and fourth-order artificial viscosity coefficients in nonlinear coefficient 
model. 

Parameters determining type of time differencing used. 

Constants in nonlinear coefficient artificial viscosity model. 

Effective second coefficient of viscosity. 

Laminar and turbulent second coefficient of viscosity. 

Effective viscosity coefficient. 

Laminar and turbulent viscosity coefficient. 

Laminar kinematic viscosity. 

Computational coordinate directions. 

Static density. 

Pressure gradient scaling parameter in nonlinear coefficient artificial viscosity 
model. 
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Symbol 


Definition 


t*x, t xyi etc. 

* 

SUBSCRIPTS 

Subscript 

j , k 

r 

t 

x,y, * 

v .£ 

T 

SUPERSCRIPTS 

Superscript 

77 

* ** 
j 


Computational time. 

Elements of shear stress tensor. 

Spectral radius in nonlinear coefficient artificial viscosity model. 


Definition 

Denotes grid location in rj, and £ directions. 

Denotes dimensional reference condition. 

Denotes differentiation with respect to physical time. 

Denotes differentiation with respect to Cartesian coordinate directions. 
Denotes differentiation with respect to computational coordinate directions. 
Denotes differentiation with respect to computational time. 


Definition 
Denotes time level. 

Denotes solution after first and second ADI sweep. 
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PROTEUS THREE-DIMENSIONAL 
NAVIER- ST OKES COMPUTER CODE - VERSION 1.0 

Volume 1 - Analysis Description 

Charles E. Towne, John R. Schwab, Trong T. Bui 

National Aeronautics and Space Administration 
Lewis Research Center 
Cleveland, Ohio 


SUMMARY 

A computer code called Proteus has been developed to solve the three-dimensional, Reynolds-averaged, 
unsteady compressible Navier- Stokes equations in strong conservation law form. The objective in this ef- 
fort has been to develop a code for aerospace propulsion applications that is easy to use and easy to modify. 
Code readability, modularity, and documentation have been emphasized. 

The governing equations are written in Cartesian coordinates and transformed into generalized 
nonorthogonal body-fitted coordinates. They are solved by marching in time using a fully-coupled 
altemating-direction-implicit solution procedure with generalized first- or second-order time differencing. 
The boundary conditions are also treated implicitly, and may be steady or unsteady. Spatially periodic 
boundary conditions are also available. All terms, including the diffusion terms, are linearized using 
second-order Taylor series expansions. Turbulence is modeled using either an algebraic or two-equation 
eddy viscosity model. 

The program contains many operating options. The thin-layer or Euler equations may be solved as 
subsets of the Navier-Stokes equations. The energy equation may be eliminated by the assumption of 
constant total enthalpy. Explicit and implicit artificial viscosity may be used to damp pre- and post-shock 
oscillations in supersonic flow and to minimize odd-even decoupling caused by central spatial differencing 
of the convective terms in high Reynolds number flow. Several time step options are available for conver- 
gence acceleration, including a locally variable time step and global time step cycling. Simple Cartesian or 
cylindrical grids may be generated internally by the program. More complex geometries require an ex- 
ternally generated computational coordinate system. 

The documentation is divided into three volumes. Volume 1, the current volume, is the Analysis De- 
scription, and presents the equations and solution procedure used in Proteus . It describes in detail the 
governing equations, the turbulence model, the linearization of the equations and boundary conditions, the 
time and space differencing formulas, the ADI solution procedure, and the artificial viscosity models. Vol- 
ume 2 is the User's Guide, and contains information needed to run the program. It describes the program's 
general features, the input and output, the procedure for setting up initial conditions, the computer resource 
requirements, the diagnostic messages that may be generated, the job control language used to run the 
program, and several test cases. Volume 3 is the Programmer's Reference, and contains detailed informa- 
tion useful when modifying the program. It describes the program structure, the Fortran variables stored 
in common blocks, and the details of each subprogram. 

A two-dimensional/axisymmetric version of Proteus code also exists, and was originally released in late 
1989. 
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1.0 INTRODUCTION 


Much of the effort in applied computational fluid dynamics consists of modifying an existing program 
for whatever geometries and flow regimes are of current interest to the researcher. Unfortunately, nearly 
all of the available non-proprietary programs were started as research projects with the emphasis on dem- 
onstrating the numerical algorithm rather than ease of use or ease of modification. The developers usually 
intend to clean up and formally document the program, but the immediate need to extend it to new ge- 
ometries and flow regimes takes precedence. 

The result is often a haphazard collection of poorly written code without any consistent structure. An 
extensively modified program may not even perform as expected under certain combinations of operating 
options. Each new user must invest considerable time and effort in attempting to understand the underlying 
structure of the program if intending to do anything more than run standard test cases with it. The user's 
subsequent modifications further obscure the program structure and therefore make it even more difficult 
for others to understand. 

The Proteus three-dimensional Navier- Stokes computer program is a user-oriented and easily-modifiable 
flow analysis program for aerospace propulsion applications. Readability, modularity, and documentation 
were primary objectives during its development. The entire program was specified, designed, and imple- 
mented in a controlled, systematic manner. Strict programming standards were enforced by immediate peer 
review of code modules; Kemighan and Plauger (1978) provided many useful ideas about consistent pro- 
gramming style. Every subroutine contains an extensive comment section describing the purpose, input 
variables, output variables, and calling sequence of the subroutine. With just three clearly-defined ex- 
ceptions, the entire program is written in ANSI standard Fortran 77 to enhance portability. A master ver- 
sion of the program is maintained and periodically updated with corrections, as well as extensions of general 
interest (e.g., turbulence models.) 

The Proteus program solves the unsteady, compressible, Reynolds-averaged Navier-Stokes equations in 
strong conservation law form. The governing equations are written in Cartesian coordinates and trans- 
formed into generalized nonorthogonal body-fitted coordinates. They are sols d by marching in time using 
a fully-coupled altemating-direction-implicit (ADI) scheme with generalized time and space differencing 
(Briley and McDonald, 1977; Beam and Warming, 1978). Turbulence is modeled using either the Baldwin 
and Lomax (1978) algebraic eddy-viscosity model or the Chien (1982) two-equation model. All terms, in- 
cluding the diffusion terms, are linearized using second-order Taylor series expansions. The boundary 
conditions are treated implicitly, and may be steady or unsteady. Spatially periodic boundary conditions 
are also available. 


The program contains many operating options. The thin-layer or Euler equations may be solved as 
subsets of the Navier-Stokes equations. The energy equation may be eliminated by the assumption of 
constant total enthalpy. Explicit and implicit artificial viscosity may be used to damp pre- and post-shock 
oscillations in supersonic flow and to minimize odd-even decoupling caused by central spatial differencing 
of the convective terms in high Reynolds number flow. Several time step options are available for conver- 
gence acceleration, including a locally variable time step and global time step cycling. Simple grids may be 
generated internally by the program; more complex geometries require external grid generation, such as that 
developed by Chen and Schwab (1988). 


The documentation is divided into three volumes. Volume 1, the current volume, is the Analysis De- 
scription, and presents the equations and solution procedure used in Proteus. It describes in detail the 
governing equations, the turbulence model, the linearization of the equations and boundary conditions, the 
time and space differencing formulas, the ADI solution procedure, and the artificial viscosity models. Vol- 
ume 2 is the User's Guide, and contains information needed to run the program. It describes the program's 
general features, the input and output, the procedure for setting up initial conditions, the computer resource 
requirements, the diagnostic messages that may be generated, the job control language used to run the 
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program, and several test cases. Volume 3 is the Programmer's Reference, and contains detailed informa- 
tion useful when modifying the program. It describes the program structure, the Fortran variables stored 
in common blocks, and the details of each subprogram, 

A two-dimensional/axisymmetric version of Proteus code also exists, and was originally released in late 
1989 (Towne, Schwab, Benson, and Suresh, 1990). 

The authors would like to acknowledge the significant contributions made by their co-workers. Tom 
Benson provided part of the original impetus for the development of Proteus , and did the original coding 
of the block tri-diagonal inversion routines. Simon Chen did the original coding of the Baldwin- Lomax 
turbulence model, and consulted in the implementation of the nonlinear coefficient artificial viscosity model. 
William Kunik developed the original code for computing the metrics of the generalized nonorthogonal grid 
transformation. Frank Molls has created a separate diagonalized version of the code. Ambady Suresh did 
the original coding for the second-order time differencing and for the nonlinear coefficient artificial viscosity 
model. These people, along with Dick Cavicchi, Julie Conley, Jason Solbeck, and Pat Zeman, have also 
run many debugging and verification cases. 
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2.0 GOVERNING EQUATIONS 


2.1 GOVERNING EQUATIONS IN CARTESIAN COORDINATES 

The basic governing equations are the three-dimensional compressible Navier-Stokes equations. These 
equations may be found in several standard references (e.g., Hughes and Gaylord, 1964; Schlichting, 1968; 
White, 1974; Anderson, Tannehill, and Pletcher, 1984.) In Cartesian coordinates, the three-dimensional 
equations can be written in strong conservation law form using vector notation as 

dQ , dE , dF , dG d ^v , dG v 

dt dx dy dz dx dy dz 


where 


Q = [p pw pv pw E r 2 


(2.2a) 


E = 


pu 

pu +p 
puv 
puw 

( [E t +p)u 


(2.2b) 


F = 


pv 

puv 

pv 2 +p 

pvw 

(E r + p)v 


(2.2c) 


G = 


pw 

puw 

pvw 

pw 2 +P 
( E r + p)w 


(2.2d) 




xy 


*Txx + VT xy+ w *x2--frr<?x 


(2.2c) 
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0 


FV = 


Re r 


‘xy 


'yy 


h yz 


(2.20 


WT jcy "h VT /y "h "^yz p r Qy 


Gv — 


Re, 


' 7 * 


UTxz + V?y Z + 


(2-2g) 


Pr, 


* 


Equation (2.1) thus represents, in order, the continuity, jc-momentum, y-momentum, z-momentum, and 
energy equations, with dependent variables p , pu, pv, pw, and £Y- 

The shear stresses and heat fluxes are given by 



In these equations, t represents time; x, y, and z represent the Cartesian coordinate directions; u, v, and 
w are the velocities in the x, y, and z directions; p, p, and T are the static density, pressure, and temperature; 
E t is the total energy per unit volume; and p, and k are the coefficient of viscosity, second coefficient of 
viscosity, and coefficient of thermal conductivity. 

All of the above equations have been nondimensionalized using appropriate normalizing conditions. 
Lengths have been nondimensionalized by L, velocities by u,, density by p„ temperature by T„ viscosity 
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by p„ thermal conductivity by k r , pressure and total energy by p r u, }, and time by L r /u r . The reference 
Reynolds and Prandtl numbers are thus defined as Re r = p,u r L r lp, and Pr r = p r u}jk,T r } 

Turbulence is modeled using the Boussinesq approach (Schlichting, 1968). The equations presented in 
this section are thus used for both laminar and turbulent flow. For turbulent flow they represent the 
Reynolds time-averaged form of the Navier-Stokes equations, with density fluctuations neglected. They 
may also be interpreted as the Favre or mass-weighted time-averaged form of the equations. With F avre 
time averaging, however, the velocities and thermal variables represent mass-averaged quantities defined by 
u = pulp, etc., where the overbar represents a conventional Reynolds time-averaged quantity. Details on 
Reynolds and Favre time-averaging procedures may be found in Cebeci and Smith (1974), and in Anderson, 
Tannehill, and Pletcher (1984). In either case, p, X, and k represent effective coefficients. For example, in 
turbulent flow p = p, + p r , where p, and p, are the laminar and turbulent viscosity coefficients, and p, comes 
from some appropriate turbulence model. The models currently available in the Proteus code are the al- 
gebraic eddy viscosity model of Baldwin and Lomax (1978) and the two-equation model of Chien (1982), 
implemented as described in Section 9.0. 

2.2 EQUATION OF STATE 

In addition to the equations presented above, an equation of state is required to relate pressure to the 
dependent variables. Any appropriate equation, or even table, could be used. The equation currently built 
into the Proteus code is the equation of state for thermally perfect gases, p — pRT, where R is the gas con- 
stant. For calorically perfect gases, this can be rewritten as 

p - (, - l)[£ r - \ f.(u 7 + v 1 2 + w 2 )] (2.4) 

where y is the ratio of specific heats, c p jc v . Here the gas constant and specific heats have been 
nondimensionalized by u}jT r . 

If the flow is such that we can assume a perfect gas with constant stagnation enthalpy, the energy 
equation may be eliminated. This assumption is reasonable, for example, in inviscid regions, and in 
adiabatic wall boundary layers if the Prandtl number is near 1 (Briley and McDonald, 1977). The stag- 
nation enthalpy is defined as 

h r =c p T+±-(u 2 + v 2 + w 2 ) (2.5) 

Here the stagnation enthalpy is nondimensionalized by i£. The temperature is thus 

T = -^[/« r -^( W 2 + v 2 + w 2) ] (2.6) 


and the equation of state becomes 

p = -^— y — Pp* 7 (i / 2 + v 2 + w 2 ) J (2.7) 

This equation of state does not require the total energy £ r , and the energy equation need not be solved. 
The total energy may be computed from 

E t — ph T — p (2-8) 


2.3 GENERALIZED GRID TRANSFORMATION 

Because the governing equations in the previous section are written in Cartesian coordinates, they are 
not well suited for general geometric configurations. For most applications a body-fitted coordinate system 


1 Note that this Prandtl number does not have a physically meaningful value, but is merely defined by a combination 

of the normalizing conditions for c p , p, and k that appear when the equations are nondimensionalized. 
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is desired. This greatly simplifies the application of boundary conditions and the bookkeeping in the nu- 
merical method used to solve the equations. The following generalized grid transformation, which can be 
orthogonal or nonorthogonal, is therefore used to transform the governing equations from physical 
(x,y, z, t) coordinates to computational (t, rj, t, t) coordinates. 

t-t(x,y,z, t) 

= *i( x ’ y> z > 0 /j q\ 

t-t(x,y,z, 0 

T = t 

In Proteus, the spatial computational domain is a cube, with t, »j, and t each running from 0 to 1. Using 
the chain rule for partial differentiation, the derivatives in the Cartesian form of the governing equations can 
be replaced using the following expressions. 


d e d_ d r d 

dx ~* x d£ + ^ dr, + dt 


_a 

dy 

_a 

dz 

_a 

dt 


a.. '’V 


— =t- 


_a_ 

dt 

d 


+ ^ dt, y dr 


d t 


JL. 

a, *=t 


, J_, r J_ 

dt Vz dr, ^ aC 

a , a r a , 
aT + ,Jf a^ + c 'aT H 


ar 


( 2 . 10 ) 


In the above equations, and in those to follow, subscripts x, y, and z, or t, r,, and t, denote partial differ- 
entiation in that coordinate direction. The only task remaining, then, is to develop expressions for the 
metric coefficients £ x , r, x , etc. In differential form we can wnte 

dt, = t x dx + tydy + t 2 dz + t,dt 
dt, = r, x dx + rjydy + r, 2 dz 4- rj/Jt 
dt = t x dx + tydy + t z dz + t,dt 
dz = dt 

In matrix form this becomes 


dt 


t x ty t 2 tr 

dx 

dr, 


*lx *ly *I Z ’ll 

dy 

dt 


t x ty t Z tr 

dz 

dz 


0 0 0 1 

dt 


Similarly, 


'dx' 


1 

£ 

J* 

A 

* 

1 

'dt 

dy 


y* y v y z 

dr, 

dz 


z l z n z t Z z 

dt 

dt 


0 0 0 1 

dz 


Therefore, 
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After taking the inverse, 


tx 


tx 




** 

x c 


*lx 

ly 

Vz 

Vt 




y c 

T t 

C* 

ty 

Cz 

Cr 


z l 

z » 

z c 

Z T 

0 

0 

0 

1 


_0 

0 

0 

1 


■«x 

ty 


ti 


V?-^ z n 

x t z n ~ Vc 

Vc-^n 

*14" 

»?x 

Vy 

»Jz 


= y 

y&-y& 

X ? Z C ~ X C Z ? 


*24 

C* 

Zy 

Cz 

Cr 

yfr-y^s 

X^-Xfr 

x s*i“V« 

*34 

0 

0 

0 

1 


0 

0 

0 

1/* 


where 

= x M z v - y v z d + -MVc - x q z J + z M y v - vc> 

*24 = *tO>{ z £ - J^) + J> T (* c Zj - *&) + z T (x^y Q - 
*34 = **(V? -J^{Z n ) ~ v«) + “ X ^«P 


and 7 is the Jacobian of the transformation, 


gtf, g. 0 

3(xj>,z) 


^ *z 

»7jc *ly *iz 

Cx ty tz 


J ~ SxOtyZz ~ *l£y) "b ^yi’lz^x Vx^z) + ^z&z&y *ly£x) 

This can be evaluated from the known physical (x, y, z) coordinates by noting J = l//~ 1 and 


X r, X i 

y Z y* J'C 

z l z n h 

j~ 1 = *{(Vc “ JW + ^0c z ? - P* 2 ;) + x e^ Z 7 _ V?) 

The metric coefficients themselves are 

£;c = *(VC 

^ y =j(X ; Z r) -X rl Z ; ) 

^=*(Vc~ 

>7x = *(>’c z £ ~ Tf z f) 

= *(^ z c - ^Z|) 

Cx = % z r,“V^ 


r \ 8{xy,z) 

a({,n. 0 


( 2 . 11 ) 


( 2 . 12 ) 


( 2 . 13 ) 
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^ = 4Vr¥-)) 

= Ax^n - 

»jf = -%-Wr z ;i; 

tt=-x£x-y£y-zJ>2 

Unless the physical coordinates (x, j?, z) are defined analytically as functions of the computational coordi- 
nates (f, r\, £), the metric coefficients must be computed numerically. The method used to do this is im- 
portant, and is discussed in Section 2.5. 

2.4 GOVERNING EQUATIONS IN COMPUTATIONAL COORDINATES 

Applying the generalized grid transformation of the previous section to equation (2.1) yields 


Qt + / + QrjVt + Q + E { f x + E ^jc + E cCjc + + + + 

" E^£ x ~~ - F^ij^ — F^Cy ~ — G v rj 2 — G^ = 0 (2.14) 

This equation is in chain-rule, or weakly conservative form. That is, the conservation flow variables are 
used, but the metrics appear as coefficients of the derivatives instead of inside the derivatives. Following 
Vinokur (1974), the strong conservation law form can be recovered by first dividing by the Jacobian then 
adding and subtracting like terms. For example, the 

Erf, 

term becomes 



Doing this for all the terms, and rearranging, results in 

+ Ff,, + + Q>r J f Frj x + F *i y + G*i z + Q*i t 1 I" + F£ y 4 * G£ z + Q£ r 


(*W 

-[ 


H 

•]-[• 


■H 

l-[ 


- r 


< -‘n *- 

E^5 X + F viy + Gytz I f ^vVx + F v rj y + G v rj 2 "| |" E v £ x + F v £ y + G V Q 2 


( 2 . 15 ) 


The last four terms, in braces, are called the metric invariant tenns. By using the expressions for the metric 
coefficients, given by equations (2.13), one can show that the metric invariants are identically zero. This is 
not necessarily true when derivatives are approximated by finite differences, however. This point is explored 
further in Section 2.5. With the metric invariant terms eliminated, no metrics or flow variables appear as 
coefficients, and the strong conservation law form of the governing equations has been recovered. 


Equation (2. 1 5) can be rewritten as 

£Q JE_ aF dG_ d¥ y dGy 

d-c + dt dv dC d£ dt, dZ 


(2.16) 


where 
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E = y(E^ + F^ + G^+Q£ r ) 
F = y (Erj x + F 'fi y + Grj z + <)*},) 
G = -Jr(E^ + F£ y + GC 2 +QQ 
E(/= y (Ej / Z x + Fy£y + G y£ z ) 
Fj/= y (E^ *] x + Fj/ rjy + Gy rj z ) 

G^=y (E K C x + F^C y + Gf/Q 
Using equations (2.2a) through (2.2g) these can be expanded as 

Q = y [p pu pv pw E-f] T 


E 


pul x + pvly + pw£ z + p£, 

(pu 2 +p)Z x + puvly + puw£ z + putt 
puv£ x + (pv 2 +p)$ y + pvw£ z + pv£ f 
puw£ x + pvw£y + (pw 2 + p)l z + pw£ t 
(E r + p)ut x + (E T + p)vZ y + (E t + p)w{ 2 + E T £ t 


(2.17a) 


(2.17b) 


A 

F 


purj x + pvr\ y + pwr\ z + prf t 
(pu +p)Tj x + puvrjy + puwrj z + putJi 
puvr] x + (pv 2 + p)r] y + pvw*] z + pv» j f 
puwr\ x + pvwrjy + (pw + p)r\ z + pwrj t 

(E t + p)w] x + (E t + p)vrj y + (E r + p)wr] z +E T ti t 



pu£ x + pvCy + pw£ z + p£, 

(pu 2 + p)t x + puv£ y + puw£ z + pidij 
puvt; x + (pv 2 + p)£ y + pvwC z + pv£ r 
puwC x + pvw^y + (pw 2 + p % 2 + pwC t 
(E r + p)ii, x + (E t + p)vtiy + (E r + p)w£ z + Ej-t.1 


E K- 


l 1_ 

J Re r 


0 


T xx%x ■*" T xy^y + T xz^z 
T xy% x "F T yy^y T yz^z 
r xz%x "F T yz^y ”F r zz^z 
Px^x fiy^y Pz^z 


(2.17c) 


(2.17d) 


(2.17e) 
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A r xx*1x r xy 1 7y T xz Y lz 

~ ~J~ r xy^x ^ T yy*1y T yz 7 lz 

T xz*b c + + 

Pxtfx PyVy "F fiz*?z 


(2.17f) 


J J r xx£x + T xy^ + T xz^z 
Gy — — r xy^x 7 yy^y T y£z 

T xJ*x T yz^y "h T z^>z 

Px^x Py^y + P£z 


(2. 1 7g) 


where 


P X = U *XX + VT ^ + WT X2 “7^-^ 

Py “ WT ;cy d" VT y>> WT yz p r 4y 
Pz~ ^ X 2 ~h V^yz “b ™t 22 4z 

In the viscous terms, the shear stresses and heat fluxes are defined exactly as in equations (2.3), except 
the derivatives in the Cartesian coordinate directions must be evaluated using the chain rule. For example, 

du _ du g du du r 

dx ~ d£ Cjc dtj Vx 3£ 

Note that F and G have exactly the same form as E, but with l replaced by >j and £, respectively. Similarly, 

A * 

i v and Gy have exactly the same form as E y, but with £, replaced by rj and C, respectively. 

2.5 METRIC INVARIANTS 

The governing differential equation in computational coordinates, equation (2.16), can be rewritten as 
/ Q \ r E? x + F^ + a 2 + QJ, ] r E^ + F^ + G^ + Qn, 1 r EC, + EC, + gc 2 + QC, -| 


(n+[ 

-[ 


EyZ x + ?vly + G vlz 
J 


r j ^ I" + F tf y + Grj z 4- Qy t j ^ ^ E£ x + F£ y + G£ g + QC, j 
J J“ Eyf} x + Tyfly + ^VVz j ^ + ? V*i>y + G V ^ 


When this equation is applied to uniform flow, E, F, G, etc., are all constant, resulting in 

i-Q, + (I)g + (E- E ^) !t ( F -r^i)^c-c,)(i) i + Q(i) { 

+ < e -^ : t)/< f - f 4t), +<g - g »' ) (t), +q (t), 

+ (E-E„)(.^-) +(F-F^^j +(0-G v )(^f^ +( *(t) =0 
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Collecting terms, 


7^ + q [(t) 


+ (E — Ey) 


+ (F~ F y) 


+ (G-gJ 



For Q t to be zero, which it should be for uniform flow, the terms in brackets must vanish. These terms 
are the metric invariants discussed briefly in Section 2.4. By using the expressions for the metric coefficients 
given by equations (2.13), one can show that in differential form these metric invariants are indeed iden- 
tically zero. When finite differences are used to approximate derivatives, however, this is not necessarily 
true. In two dimensions, when the centered difference formula of equation (5.1) is used to approximate 
derivatives, the metric invariants do turn out to be identically zero. But in three dimensions, when the 
metric coefficients are computed numerically using equations (2.13), the metric invariants are not identically 
zero. 

To show this, let M 2 denote the second metric invariant term (the second bracketed term in equation 
(2.18).) Then, using equations (2.13) for the metric coefficients and applying difference operators, 

m 2 ~ <^( v 6 c? ~ V V) + W V ~ V s t z ) + W V ” V 

Without loss of generality we can let A£ = Arj = AC = 1/2. Then, using central differences, 

M 2 — <^[0 ; y + l ~ yj - l )( z k + 1 ~ z k - 0 ~ + \~ ^k- 0( z j + 1 z j — l)-l 

+ + i ~yk- l)( z f + l ~ z i- l) ” 0>/ + i — J'/- l)( z £ + 1 z k - l)] 

+ -j- 1 ~ yi - 1)(^* +i ~ z j - 1) ~ (ty* + 1 ~ 1 X 2 / + i z /- 1)] 

The subscripts /, y , and A: represent grid point indices in the £, and £ directions. For notational conven- 
ience, terms without an explicitly written z,y, or A: subscript are understood to be at z,y, or Expanding, 

M 2 — + i 7 * + 1 ” J'y - i 2 * + l + l 7 * - 1 yj - I 7 * - 1) “ + \ 2 j + 1 — ¥k - l 7 y + 1 — ¥k + i 7 / - l - i 7 y - i)H 

+ <5, E(y* + i z t + 1 — y* - i^t + 1 “.y* + i 7 * - 1 + yk - & - 1) ~ + i 7 * + 1 — ^ - i 7 * + i — ^ - i 7 * - ^ 

+ ^C(X|* + I^+ I — >V - l 7 / + 1 l 7 ;- 1 - l 7 /— l) — 0y + l Z 4 + 1 1 2 I + 1 + 1 7 <- 1 

Performing the final difference operation, 

M2 ^ (y,+ i* k +i —yj-i z k+\ ~yj+\ z k-\ +>/- \ z k- j)i+ 1 “ (yj+i 2 k+ 1 — .b - 1 7 *+ 1 “^y + 1 7 *- 1 + ^y - i 7 * - iX - 1 
~ (V* + i 7 y + 1 — - i 7 y + 1 — yk + \ z j - 1 +yk- \ z j- lX + 1 + O'* + i 7 y + 1 ” - l 7 y + 1 + l 7 ; - 1 ^ - 1 2 j - A - 1 

+ (y*+ l 7 * + 1 — l 7 i + 1 + 1 7 4 - 1 +yk- \ Z i- l)y+ 1 — + 1 7 I + 1 ~ 1 7 i + 1 1 7 4 - l + — l 7 i - l)y - 1 

~ (y, + + 1 — yi - \ Z k + 1 “ ^ + 1 7 A - 1 + ^4 - 1 7 A - l)y + 1 fri + 1 7 A + 1 “ - 1 7 A + 1 “ ^4 + 1 7 A - 1 X - l 7 * - lV - 1 

+ (y r + \ 2 j + j _ + 1 - > t - + \ z j _ i +^ _ i^,- _ i)* + 1 — + i 7 y + 1 “ - i z j + 1 “ -^4 + i z j - 1 + - i 7 y - i>A - 1 

Finally, collecting terms. 
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A/ 2 — J>/ + I J + \Jc( z i+ l,jjc+ 1 z i+\,jJc-\ z i,j + \Jc + \ + z i,j + \ Jl — V 

+ K+ \ J- ljc(~ z i+ \JJc + 1 + z i + IJJc - I + z i,J- z i,j-lJc-\) 

+ - 

For M 2 to be identically zero, the terms in parentheses in the above equation would have to vanish 
identically. This is clearly not the case for a general three-dimensional coordinate system. 

There are fixes that have been developed to ensure that the finite difference equations do satisfy uniform 
flow, including the use of averaging formulas for the metric coefficients or simply subtracting the error from 
the equations (Pulliam and Steger, 1978). These methods are somewhat inelegant, however, and can be 
expensive to use. A cleaner and completely rigorous method is to rewrite the formulas defining the metric 
coefficients, equations (2.13), in conservation form (Thomas and Lombard, 1979). Using s« as an example, 

£W0y*-J W 

= J l(y n z \ - O'c 2 ),] 

In Proteus , therefore, the metric coefficients are actually computed using the following equations. 

~ ^ 2 )r,] 

Zy = Jlix&r, - 

- OyOr,] 

Vx — — (V{ 2 ){] 

Vy = ~ (*{ 2 )jl 

n z = JU. x t?)i ~ (*{?){] ( 2 - 1 9 ) 

t: x = Jl(y>z) r} -(y n z) l ] 

Z 2 = JUx^n - ( x ^l 

Zt=-x^ x -y^y- z £z 

r I:~~ X z r lx y^y ™“ 

C/ = - x £ x -y£ y - z t C 7 

To verify that computing the metrics in this way does lead to metric invariants that are identically zero, 
we now reevaluate M 2 , the second bracketed term of equation (2.18). Applying difference operators, 

M 2 si S { [<5£ (y v z) - S' n (y ? z)] + 6 n [8\ tyz) - <5£ (y^z)] + 8^S' V (y^z) - <5£ (y^z)] 

Note that a distinction is made, for now, between the difference operators outside the brackets, S ( , <5,, and 
<5 { , and the difference operators inside the brackets, S', , <5' , and S\ . The operators outside the brackets 
represent derivatives of the metric coefficients in M 2 . These terms originated as part of the flux terms in 
equation (2.16). The operators S { , <5,, and S ( are therefore the same as the operators used to represent de- 
rivatives in the governing differential equation. Now, note that the finite difference approximation of M 2 
will vanish identically if, for example, 


This will be true if S', = 6 ( , S', = <5„ and S[ = <5 ( . This can be verified by expanding <5 { (<5 { /) and <5 { (<5 { /) 
using the centered difference formula of equation (5-1), and comparing. When the metric coefficients are 
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computed using equations (2. 19), therefore, the derivatives of the parenthetical terms must be approximated 
using the same difference operators as those used to represent derivatives in equation (2.16). It does not 
matter how the x 0 x,, etc., inside the parentheses are computed. 

This procedure for computing the metrics ensures that the last three metric invariant terms in equation 
(2.18) are identically zero when differenced. The first metric invariant term must be handled somewhat 
differently. Setting it equal to zero gives 

(I)/(t) ! + (t)/(t) ! =» < 2 - 20) 

This is a statement of the geometric conservation law described by Thomas and Lombard (1979). For grids 
that do not change with time, this equation is, of course, automatically satisfied when differenced. However, 
for time-dependent grids it is not. In that case, the grid transformation Jacobian J should be found by 
solving equation (2.20) at the new time level, using the same differencing scheme as in the governing flow 
equations, and not computed algebraically from equation (2.12). The current version of Proteus does not 
solve this equation, and thus strictly applies only to time-independent grids. 
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3.0 TIME DIFFERENCING 


The governing equations are solved by marching in tune from some known set of initial conditions using 
a finite difference technique. The time differencing scheme currently used in Proteus is the generalized 
scheme of Beam and Warming (1978). The time derivative term in equation (2.16) is written as 


8Q AQ” _ £i d(AQ") 1 dQ n 0 2 AQ n 1 

5t _ At 1 + 0 2 5t 1 + ^2 1+^2 At 

+ o [( 01_ T“ 02 ) At,(At)2 ] 


or, 


A „ 

AQ” = 




d(AQ n ) 


1 + 9~) 


dr 


+ 


At 


dQ n 


l+0o dz 


+ 


1 + 0, 


aq” - 1 + oj 


[( 0 ,-- 0 2 )(Av) 2 , (At) 3 J 


(3.1) 


where AQ' 1 = + 1 - Q". The superscripts n and n -h 1 denote the known and unknown time levels, re- 

spectively. 

The parameters 6 X and d 2 deteimine the type of time differencing scheme used. Some of the methods 
available with the above formula are given in the following table. 


*1 

^2 

Method 

Truncation Error 

0 

0 

Euler explicit 

O(At) 2 

0 

-1/2 

Leapfrog explicit 

O(At) 3 

1 

0 

Euler implicit 

O(At) 2 

1/2 

0 

Trapezoidal implicit 

~‘(At) 3 

1 

1/2 

3-point backward implicit 

O(At) 3 


Note that even though the generalized time differencing formula includes explicit methods, the Proteus code 
assumes an implicit method is being used. Note also that the truncation error listed in the table is the error 

in the expression for AQ". The overall numerical method used in modelling the differential equations re- 
quires AQ'VAt, so the order of the overall method is this truncation error divided by At. 

Solving equation (2.16) for dQjdr and substituting the result into equation (3.1) for d(AQ n )jdr and 
dQ'Ydr yields 

a„ 0jAt f d(AE n ) <3(AF*) d(AG n ) \ Ar f dE n dE* dG n \ 

AQ ” ~ 1 + 0 2 ( dt + dr, + 0C J 1 + 0 2 \ di + dr, dc J 

0! At f d(AE y n ) d(AF v n ) d(AG y n ) \ At / dG v n \ 

+ i + 0 2 ( at + + di : J + i + 0 2 ^ at. at, * J 

+ -p^r- aq” - 1 + oT(«, - \ - e 2 )(it) 2 , (*r) 3 J (3.2) 
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4.0 LINEARIZATION PROCEDURE 


4.1 INVISCID TERMS 

A A A A 

Equation (3.2) is nonlinear, since, for example, AE rt = E n+1 — E* and the unknown E n + I is a nonlinear 
function of the dependent variables and of the metric coefficients resulting from the generalized grid trans- 
formation. The equations must therefore be linearized to be solved by the finite difference procedure used 
in Proteus . This is done by expanding each nonlinear expression in a Taylor series in time about the known 
time level n. Letting G represent any nonlinear expression, 

G n + 1 = G n + ( ) At + O(Ar) 2 (4.1) 


where 

dG dG_ dp dG g (pu) dG d(pv) gg g (pw) dG dE T 

dr ~ dp dr d(pu) dr d(pv) . dr d(pw) dr dE r dr 


Note that for linearization purposes only the metric scale coefficients have been assumed to be locally inde- 
pendent of time. Note also that for this linearization procedure to be second order accurate, dGjdr (and 
therefore dp I dr, d{pu)jdr, etc.) need only be first order accurate. Using forward differences, then, so that 


(*)- 


n n + 1 n 
P -P 


At 


+ 0 (At) 


A P 
At 


+ 0(At) 


etc., equation (4.1) becomes 


^ n + 1 


G n + 




( % + ( iw y^f + ( )" 4(,>v) " 


(4.2) 


As an example the d(puv£ y )ld£ term from the x-momentum equation (part of the second element of 

3E/3£) will be used. The nonlinear part of this term is (pwv)” + '. Rewriting this in terms of the dependent 
variables, 


(puv) 


n + 1 


(p“)(p v ) 

P 


Using equation (4.2), this is linearized as 

(puv) n + 1 = (puv)” - (uv)”(p” + 1 - p”) 4- v”[(pu)” + 1 - (pu)"] + u”[(pv)” + 1 - (pv)”] + O(At) 2 
which can be rewritten as 


A(puv)” = — (uv)”Ap” + v”A (pu)™ + u n A(pv)” + 0( At) 
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A 


This linearization procedure, when applied to the entire AE n term in the vector equation (3.2), can be 
written as 


AE 


n 


) n 

AQ n + 0( At ) 2 


(4.3) 


A A 

where (<3E/<3Q) n is a Jacobian coefficient matrix (not to be confused with the Jacobian J of the generalized 

A A 

grid transformation.) Similar equations can be written for AF n and AG*. 

AAA 

Each term in each element of E, F, and G, given by equations (2.17b) through (2.17d), is linearized using 

A A A A 

the above procedure to generate the elements of the Jacobian coefficient matrices 5E/5Q, dFjdQ, and 
dGjdQ. (Note that 5E/5Q = Jd E/3Q.) When this is done <5E/3Q can be written as 


dE 

dQ 


It 


dp 

dp 


dp 
dp 

ZR-* _ 

dp 


ly “ V A 




I* 




dp 

d(pu) 


lx 


- 4 -*) 


dp 

v£x + ~d<JU) ** 
dp 

w£x + £x 

dp 

J**+*Tk 


dp 

" £ > + ~d{^) £x 
dp 

+ 3/ -.A 


5, 


" £x + -afe- £x 

9p 

v 5z + -57rzrf. 




Jily +/l 


ao?v) 

dp 

d(pv) 


It +/l + + 

/A+/1 


dp 


d(/>w) 

dp 




dp 

d£ r 

dp 

d£ r 

dp 

a£ r 


lx 

i, 

S. 


d{pw) 


£ ' +/l ( 1 + 


( 4 . 4 ) 

A A A A 

where fi = u£ x + v£ y + and f 2 = (E T + p)jp . The Jacobian matrices 5F/3Q and dGjdQ have the same 

A A 

form as dFjdQ, but with £ replaced by rj and £, respectively. 


The linearized pressure terms have deliberately been left in terms of dpi dp, dpjd(pu), etc. The ex- 
pressions to be used for these derivatives depend on the equation of state. Those currently built into the 
Proteus code, for a perfect gas, are presented in Section 4.3. 


4.2 VISCOUS TERMS 


AAA 

The nonlinear viscous terms in equation (3.2), involving AEk, AFv, and AGk, must also be linearized. 

A A A 

To do this, the elements of Ek, and G^, given in equations (2.17e) through (2.17g), must first be re- 
written in terms of the dependent variables, and with derivatives in the Cartesian directions transformed to 
derivatives in the computational directions using the chain rule. When the resulting expressions are sub- 
stituted into equation (3.2), mixed second derivatives appear as well as second derivatives in a single coor- 
dinate direction. The mixed, or cross, derivative terms would lead to considerable complications in the 
implicit numerical solution algorithm if they were linearized using the procedure presented in Section 4. 1 . 

A A A 

The two types of second derivatives are thus treated differently, and E^» Fy y and Gv are written as 


E y — Ej/^ 4* 

Fy=F Vi + Fy 2 (4.5) 

AAA 

G^= + Gy 2 
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where t Vv ¥ Vv and G Kl only contain derivatives in the rj , and C directions, respectively, and E K2 , Fk 2 , and 

Gv 2 contain derivatives in the remaining directions. The fully expanded expressions for E Vli Ek 2 , etc., are 
fairly long, and therefore are presented in Appendix A. 

4.2.1 Non-Cross Derivatives 

Examination of the elements of Evi in equations (A. 2a) through (A. 2d), and (A.2f), shows that eveiy 
term has the form fg^ where g is a function of the dependent variables, and /is a function of p, 2, k , and/or 
the metric coefficients. Expanding in a Taylor series about time level n gives 

(/*{)*' + 1 = (/*?)" + [- ] At + 0(Ar) 2 

For linearization purposes only, we will assume /is locally independent of time. We can thus write 
where 

dg dg dp dg d(pu) [ 

dr dp dr d(pu) dr 

Therefore 


(fg^ +l =(fg^+r-^ 


d g . . 3g 

dp Ap + d(pn) 




A(pw) + — + O(At) 


As with the inviscid terms, the linearization procedure for the entire AE y " viscous term in equation (3.2) 
can be written as 


AE n u = 


8Ey x 

A 

dQ 


AQ” + O(At) 2 


(4.6) 


A A 


Similar equations may be written for AF^j and AGki- The Jacobian coefficient matrix dE^/dQ is 
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where 

a xx = (2^ + 

+ ( 2 M + ^){y 2 + M^z 2 

a zz = M^x 2 + /*$/ + (2 m + k)Z? 
a xy =((i + W x Z y 
a x2 ~ (M + Wxtz 
a y2 = (n + W y Z z 

«0 = *«x 2 + e/ + ^ 2 ) 
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Like the pressure terms discussed earlier, the form of the temperature terms will depend on the equation 
of state being used. Those currently built into the Proteus code, for a perfect gas, are presented in Section 


4.3. 


Note that in equation (4.6) the derivatives appearing in the Jacobian coefficient matrix dEyJdQ are also 
to be applied to the AQ" appearing outside the parentheses. For example, the element in the second row 

and second column of 3E vJdQ, which corresponds to the A (pit) term in the x-momentum equation, is 
a xx d(\ lp)Id£. For this term, the notation used in equation (4.6) means 

n 
22 

n d ( A (pufjf \ 

axx aM p" I 



The Jacobian coefficient matrices for the remaining non-cross derivative viscous terms, dFyJdQ and 
dGyJdQ , have the same form as 5E vJdQ, but with £ replaced by r\ and £, respectively. 

4.2.2 Cross Derivatives 


As stated earlier, linearizing the cross derivative viscous terms in the same way as the remaining terms 
is very complicated within the framework of the implicit numerical solution algorithm used in Proteus . 
They are therefore simply lagged (i.e., evaluated at the known time level n and treated as source terms.) 
As noted by Beam and Warming (1978), this does not lead to a formal accuracy loss since 

AE \ = AEy~ 1 + 0( At ) 2 

AF^ = AF£~ 1 + 0( At) 2 (4.8) 

AGy 2 = AG n y ~ 1 + O(At) 2 
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4.3 EQUATION OF STATE 


The expressions to be used for dp/dp , dTjdp, etc., which arise from the linearization procedure, depend 
on the equation of state. The equation currently built into Proteus is for perfect gases, and can be written 
as 


or, in terms of temperature, as 


p = (y- l)[ £ 7 '-y p(“ 2 + v 2 + w 2 )] 

T “i[-£-T(“ 2+ ' a+ " 2 )] 


f UVl 1 v Ci Vi V 


With this equation of state, then, the appropriate 

dP y —1 f 2 , 2 , 2 "\ 

ip — — +v +w ) 


dT 

dp 


dp 

(y-i)w 

d{pu) 


dp 

(v- i)v 

d{ P v) 

dp 

(y - l)w 

d(pw) 


Cl> _ 

II 

v- 1 

1 

’ Et 

1 / 2 , : 
~p( u + v 

Cv 

2 

L p 


dT 

u 


d(pu) 

CyP 


dT 

V 


d(pv) 

CyP 


dT 

w 

d(pw) 

CyP 


dT 

1 


ll 

tq 

fO 

CyP 


(4.9) 

(4.10) 

(4.11a) 
(4.11b) 
(4.11c) 
(4. lid) 
(4. lie) 
(4.12a) 

(4.12b) 
(4.12c) 
(4.1 2d) 
(4.12e) 


If constant stagnation enthalpy is assumed, as discussed in Section 2.2, the appropriate equation of state 
is 

p= ~ y 1 (M 2 + v 2 + w 2 )j (4.13) 

and the temperature becomes 

7 ' = ^[ /l7 '-y ( « 2 + v 2 + w 2 )] (4.14) 

With these equations, the derivatives of p and T with respect to the dependent variables are 
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(4.15a) 


dp_ 

dp 


| ~h T + y (« 2 + v 2 + w 1 ) ] 


Sp _ y - 1 

d{pu) - y u 

dp y — 1 

~d(ri = ~ v 

dp y - 1 

3(pw) ~ V W 


37 

dp 


1 

CpP 


(u 2 + v 2 + w 2 ) 


37 u_ 

3(p«) 

37 v 

3(pv) c pP 


37 w 

3(pw) c pP 


(4.15b) 
(4.15c) 
(4. 1 5d) 
(4.16a) 
(4.16b) 
(4.16c) 
(4.16d) 


4.4 LINEARIZED GOVERNING EQUATION 

The linearized form of equation (3.2) can now be written as 



At / 3E | 3F ( 3G j | At 


3El 


A A 

dFv dGi 


1 +0 2 \ d£ dr, dZ 


1 + 3 , 


■ + 


+ 


(1 +3 3 )At 
1 + $2 

e 2 


3Et/ 3F l 


dt dr, 


d G v 

~dT 


d£ dr, dc 


3 3 At ( dEy 2 v 2 + ^v 2 


n — 1 


1 + 0 2 \ dl dr, 


dC 


1 + 


AQ n “ 1 + O (e, - i- - 9 (At) 2 , (9 3 - 0 ] )( At) 2 , (At) 3 


(4.17) 


There are a couple of things that should be mentioned about this equation. First, this equation is in 

A A 

so-called "delta" form. We will actually be solving this equation for A Q" and recovering Q"* 1 from 

Q* + 1 = A Q" + Q\ And second, in the coefficients of the cross derivative viscous terms the time differencing 
parameter 9 1 has been replaced by 0 3 . For second order time differencing (i.e., if 0i = d 2 + 1/2), 0 3 should 
be set equal to 6\, For first order time differencing, however, 0 3 can be set equal to zero without losing 
accuracy. 


Proteus 3-D Analysis Description 


4.0 Linearization 31 


PAGE. 


.WTEMOONALtf BLANK 



5.0 SPACE DIFFERENCING 


To solve equation (4.17) an evenly spaced grid is defined in the computational (£, r/, Q coordinate sys- 
tem. Spatial derivatives are then approximated by finite difference formulas. First derivatives in the £ di- 
rection are approximated using the following second-order central difference formula. 


/ df \ -£+i,y,fc fi-\J,k 

[ dt ) - VU k - 2A£ 


(5.1) 


The subscripts i,j, and k represent grid point indices in the s , r\, and C directions. The computational grid 
spacing A£ is constant, and equal to l/(JVi - 1), where N t is the number of grid points in the £ direction. 
Similar formulas are used for first derivatives in the >/ and C directions. 


The non-cross derivative viscous terms in the c, direction in equation (4.17) all have the form 
where Q represents one of the elements of Q. Using central differences this is approximated by 


-&\f-§r (^0] * * 

L — h',y, k 

""AT k 6 & A ®‘+ 1 / 2 ,y , k fi— 1 / 2 ./, 1 / 2 , y, ft) 

= — w+ , temj.k] 

- A _ 1/2J , *[(yAfi) /Jt *- (gA0 ( -_ *]} 

= {Of,/, * + 7l + 1 y fc)[(gA0, + 1 ,y, * _ (?^0i,y, 

2(A?r 

— k 7/- l,y, £ — — 1 J, /J) 

= T — 1 k fij, k)(S^Q)i — 1 , 7 , £ 

2(Af) 

— (7/ — l , 7 , £ 2 fij ; £ + — 1 ,y, £ 

+ ( Ay, £ + f + 1 J, +1,7, 

Similar formulas are used for second derivatives in the r\ and £ directions. 

Cross derivative viscous terms in the %-r] direction are evaluated using the following central difference 
formula. 
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A /}* 


- S i(f 6 r,S)iJ,k 

U, k 

— '2A~i + hV» + 1 J,k~ fi— \J, ki^tj &)i — 1 ,/, 

= 4A £&ri ^ + W> + W + l>k~8i+ l,y — 1, A:) 

— — 1 ,y, — \J+ \ , k ~" &i — 1 ,y — I » (^-^) 

Similar formulas are used for the remaining cross derivatives. Note that this formula is only needed for the 
source terms, since the viscous cross derivative terms are lagged one time level. 


When first derivatives are needed normal to a computational boundary, such as for Neumann boundary 
conditions, either first- or second-order one-sided differencing is used. The first-order formula at the f = 0 
boundary is 




1 J.k 


- (flj, k~ f\ ,y, k) 


(5.4) 


and at the £ = 1 boundary, 



) - (AW. k~ f 1 J, k ) 

' N x ,],k 


The second-order formula at the £ = 0 boundary is 



1 

2A£ 


( ^f\ J> k + k fl,j,k) 


and at the f = 1 boundary, 

2A^~ ^ N \ 4 /v, - 1 J, k + 3 A, y, ' 

Similar formulas are used at the computational boundaries in the r\ and £ direction. 



(5.5) 


(5.6) 


(5.7) 
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Choosing boundary conditions is perhaps the most important step in solving a flow problem with 
Proteus. Since the equations being solved at interior points are the same for every problem, the boundary 
conditions are what determines the final flow field for steady flows. 

With the difference formulas presented in Section 5.0, N eq boundary conditions are required at each 
computational boundary, where N eq is the number of equations being solved. Note, however, that this is 
a numerical requirement, not a mathematical one. For example, for one-dimensional Euler flow N eq = 3. 
However, characteristic theory shows that, mathematically, only two conditions may be specified at a sub- 
sonic inflow boundary, and only one at a subsonic outflow boundary (Pulliam, 1986a). Some sort of ex- 
trapolation is typically used for the additional numerical boundary conditions. 

A variety of boundary conditions are built into the Proteus code, including: (1) specified values and/or 
gradients of Cartesian velocities u , v, and w, normal velocity V n , coordinate direction velocities V 0 V q , and 
Fj, pressure p, temperature T, and density p; (2) specified values of total pressure pr , total temperature 7Y, 
and flow angles; (3) linear extrapolation; and (4) spatial periodicity. Another useful boundary condition is 
a "no change from initial condition" option for u, v, w, p, T, p, pr, and/or 7Y- Provision is also made for 
user- written boundary' conditions. The boundary conditions may be steady, unsteady, or time-periodic. 
The exact combination of boundary conditions to be used will depend on the problem being run. 

The boundary conditions in Proteus are treated implicitly. They may be viewed simply as additional 
equations to be solved by the ADI solution algorithm. And, in general, they involve nonlinear functions 
of the dependent variables. They must therefore be linearized using the procedure described in Section 4.0. 
The following sections describe this linearization for the general types of boundary conditions currently built 
into Proteus . 

6.1 NO CHANGE FROM INITIAL CONDITIONS, Ag=0 

This boundary condition simply sets the boundary value of the function - equal to its initial condition 
value. It can be written as 


A/ = / + 1 -V = 0 (6.1) 

A 

In general, g can be a nonlinear combination of the dependent variables Q. Linearizing g using the proce- 
dure described in Section 4.0, we get 


g 


n + 1 


= g n + [jl-'\ AQ” + O(Ar) 7 


(6.2) 


Neglecting the 0( At ) 2 linearization error, the linearized form of equation (6.1) can thus be written as 


dg 

A 

dQ 


AQ n = 0 


(6.3) 


6.2 SPECIFIED FUNCTION. g=f 

A specified function at a boundary can be written simply as 

g n + ] =f (6-4) 
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where g is the function being specified and / is the value being specified. Note that / can vary along the 
boundary, and can be time-dependent. Using equation (6.2) and neglecting the linearization error, the 
linearized boundary condition becomes 


AQ n =f-g n (6.5) 


6.3 SPECIFIED COORDINATE DIRECTION GRADIENT, deld<t> = f 


A specified gradient of a function in a coordinate direction can be written as 


d<j> 1 J 


(6.6) 


where g is the function whose gradient is being specified, f is the specified value, and <j> is the coordinate 
direction £ , »/, or £■ Note that /can vary along the boundary, and can be time-dependent. 

The linearized form of g is given by equation (6.2). The linearized form of equation (6.6) can thus be 
written as 


dg 

d<j> 




=/+ 0 ( At) 2 


(6.7) 


Replacing differential operators with difference operators and neglecting the linearization error, the 
linearized boundary condition can be written as 


S 


<t> 



=/- v* 


(6.8) 


..here &$ represents the one-sided difference operator to be used at the boundary. Options are available in 
Proteus to use either first-order two -point or second-order three-point differencing. 

Note that this boundary condition is a specified value of the derivative with respect to the computational 
coordinate, not with respect to the physical distance in the direction of the computational coordinate. 
Following Kom and Korn (1968), and using the properties of the generalized coordinate transformation, it 
can be shown that for the £ direction the two derivatives are related by 

dg 1 gg 


where 

% = 4- [(>!/, - + fa/x - nlf + (v£y - n ytx) 2 ] 


Similarly, for the r\ direction, 


dg 

ds i 


1 dg 


where 
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G m - jr W- - «/ + «A>- + <«, - Vx) 2 ] 


And, for the £ direction, 





dC 


where 

(j^r = -J2 ~ ?z^>) + (^z f !x ~~ Zx*lz) + i^x*1y — €y*1x) ] 


If the value /= 0, of course, the two derivatives are equivalent. 

6.4 SPECIFIED NORMAL DIRECTION GRADIENT, Vg. n =/ 

A specified gradient of a function normal to the boundary can be written as 

Vg n + 1 . „ =/ (6.9) 

where g is the function whose gradient is being specified, /is the specified value, and n represents the unit 
vector normal to the boundary. Note that / can vary along the boundary, and can be time-dependent. 

For illustrative purposes, assume we are specifying a gradient normal to a constant £ boundary. Then 

n ~ | | — ~m i T- £yj + m 

where 

= V^JC 2 + + ^2 


Equation (6.9) can then be written as 


1 / n + 1 x. i + 1 t i — ^ 1 t \ f 

-TFT fee %x 8y + 


(6.10) 


Using the chain rule to expand + \ and + 


n + 1 n+U , « + I , * + V 

fo ^x *lx “b ££ ’X 

gy ~g\ iy+gq Vy+gQ Cy 

n + 1 n+ U , « + 1 i_ cr” V 

& = £f <=2 + ^ ^2 + ^C ^ 


Substituting into equation (6.10) and rearranging, 

g\ + '(^ 2 + + ^ 2 ) + + '(^x + + Zztz) + g^ (%x£x + %y^y + %£z) = m f 


Solving for g? + 

(|- )” + ' - 4- &-X + + 5*)(|- 


+ 1 



(^x+ZyZy + tztz) 



( 6 . 11 ) 


Now, in order to incorporate this equation into the ADI solution procedure used in Proteus , the dgjdrj and 
dgjdi terms in equation (6.1 1) are lagged one level, and evaluated at time level n instead of n + 1. Strictly 
speaking, this introduces an 0 (At) error into the solution. In practice, however, the actual error will depend 
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on the degree of nonorthogonality of the coordinates near the boundary. For orthogonal coordinates no 
error is introduced. 


Using equation (6.2), and introducing difference operators and neglecting the linearization error, we can 
now write the linearized boundary condition as 


dg 

A 

dQ 


AQ 


m 


m 


«x*»x + t/ly + V" - “T + + ( 6 - 12a ) 

m 


where S < represents the one-sided difference operator to be used at the boundary. Options are available in 
Proteus to use either first -order two -point or second-order three-point differencing. 

Note that the unit vector n in equation (6.9) is in the direction of increasing £. Therefore, a positive 
value for/in equation (6.12a) indicates a flux in the direction of increasing £. Thus, a positive/at £ = 0 
implies a flux into the computational domain, and a positive / at £ = 1 implies a flux out of the computa- 
tional domain. 


Specifying a gradient normal to a constant y\ or £ boundary is done in an exactly analogous manner. 
The resulting equation for an rj boundary is 


dg 

A 

dQ 


AQ" 


= 7^f V (rixtx + Vyiy + 'Iz^z^^g” + *ly£y + - ( 6 - 1 2b ) 

m m 


where 


m = y/rjx 2 + r\y + r\z 


For a C boundary the equation is 



4- V + Cjr^ + «*)*«*" T Kjflx + tyly + C A) V” “ V" ^ 6 - 12c ) 

at? 


where 

A positive value for /in equation (6.12b) indicates a flux in the direction of increasing n* Thus, a positive 
/ at >/ = 0 implies a flux into the computational domain, and a positive /at rj = 1 implies a flux out of the 
computational domain. 


6.5 LINEAR EXTRAPOLATION 


Linear extrapolation from the two adjacent interior points is also available as a boundary condition. 
At the £ = 0 boundary, where / = 1, this can be written as 


Si 


+ 1 


T ■ 


(6.13) 


Note that this is equivalent to setting (d 2 gld£ 2 );+i = 0. Using equation (6.2), we can wnte the linearized 
boundary condition as 


Sg 

A 

dQ 


AQf - 2 


dg_ 

A 

dQ 


AQ"+ , + 


i + 1 


dg 

A 

dQ 


aqT + 2 — tf + 2tf+i-af+2 (6 - 14) 


/ + 2 


38 6.0 Boundary Conditions 


Proteus 3-D Analysis Description 



Analogous extrapolation boundary conditions can easily be written for the remaining boundaries. 
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7.0 SOLUTION PROCEDURE 


7.1 ADI ALGORITHM 


The governing equations, presented in linearized matrix form as equation (4.17), are solved by an alter- 
nating direction implicit (ADI) method. The form of the ADI splitting is the same as used by Briley and 
McDonald (1977), and by Beam and Wanning (1978). Although the split equations can be developed in 
more than one way, in this discussion the approximate factorization approach is used. 

Letting LHS(4.17) represent the left hand side of equation (4.17), we can write 



where I represents the identity matrix. Note that in this equation, using the dj dc, 
notation used is meant to imply 

J_( 3E 
« ^ 5Q 3Q 

The term in braces in equation (7.1) can be factored to give 




term as an example, the 



The last two terms represent the splitting error. Note that, since AQ" = O(A-r), these terms can be neglected 
without affecting the overall time accuracy of the algorithm, even when second order time differencing is 
used. 
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Equation (4.17) can thus be rewritten in spatially factored form, and, neglecting the temporal truncation 
and splitting error terms, becomes 


l + 


d / 3E 


dE v 


i + e 2 3M 


*|At did F 


dT v 


1+^2 dr 1 \ £Q dQ 


1 + 


d dG 


dG v 


At / 3E 3F 3G \ At 

i + * 2 \ ac y i + 1 


dEy x + dVy x + dG Vi ^ + (l + 0 3 )A 


1 + aQ ao 


3F l 1- <3Gp 


AQ 


0 2 \ *1 % 


l+0 2 \ # " d v T 3C 


0,At / 3E K 3F k 3G v 


1 + 0 2 \ d( 




( 7 . 3 ) 


Equation (7.3) can be split into the following three-sweep sequence. 
Sweep 1 (C direction) 


aq* + - 0,At d 


i + e 2 at 


dE ' AQ* 


dQ 


At d 

\+e 2 dt 


dE Vi 

A 

dQ 


A, 

AQ 


At f d E.dF dG V* At ( dE ^i dEv i dGy i 

+ -T 'r—^r I + , I — ~ 1 7 r 


i + e 2 \ di d v dc I i + e 2 \ di d v ac 


(1 + 0 3 )At f dEy 2 dFy 2 dGy 2 

+ \+e 2 \ dk + dr, + ac 


0oAt ( <3Ej/ d¥ v dGp/ 


2 +-^+- 


a „ » - 1 

r, 

2 


l + 0 2 \ d£ dtj 




do 


+ 


1 + 6-1 


AQ 


n - 1 


(7.4a) 


Sweep 2 (n direction) 


A** 

AQ + 


6jAt ^ 
1 + $2 


A 

3F 

A 

aQ 


A** 

AQ 


^iAt q 

1 + 0 2 dt. 


a Fy t 

A 

aQ 


A** 

AQ 


a « 

= AQ 


(7.4b) 


Sweep 3 (C direction) 


a 0,A? q 
AQ + • ° 


i + a 2 ac 


A 

ac 

A 

aQ 


A M 

AQ” 


0,At q 

i + 0 2 ac 


A 

aG t 


aQ 


AQ” 


£.** 
= AQ 


(7.4c) 


In the above equations, Q* and Q” represent intermediate solutions to the governing equations. 2 It should 
be noted that in Proteus , physical (i.e, n + 1 level) boundary conditions are used during the first two ADI 


2 


The notation here is somewhat inconsistent. The quantity AQ" = Q" +1 — Q", but AQ =Q — Q n , not 
Qn + i _ q\ Similarly, AQ” = Q” — Q", not Q' 1 + 1 — Q . 
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sweeps. This introduces an 0 (At) error in <3Q/dr on the boundary for unsteady flows, but no error for 
steady flows. This point is discussed in detail by Briley and McDonald (1980). 

Applying the spatial differencing formulas of Section 5.0 results in 

Sweep 1 (c direction) 


AQ; + 


0,At 


(1+0 2 )Ac 


A 

dE 

A 

dQ 


AQ i _ 1 -f(2a-l) 


/- 1 


A 

dE 

A 

dQ 


AQ, +(!-«) 


A 

dE 

A 

3Q 


A* 

AQ i+ i 


/+ 1 


0jA r 


■ [(fi - 1 + fifii- iAQ*_ 1 - (fi - 1 + 2fi+fi+\) n gi‘ A Qi + Ui' + fi+ i)V +1 AQ /+ 1] - 


(1 + 0 2 )(A^ 

— (<5 { E + + <5 C G)” + (<3|E Ki + 5 v Y Vx + <5^)" 


+ ° L + +^ ~ + <W 2 + dfi v ) n - + <5^ + 


\«-l 


1+0: 


AQ 


/? — 1 


(7.5a) 


Sweep 2 (n direction) 
0 ,At 

(TTW 


0j At 


A** 

AQ; + 


-^T ) AQj: , + (2a - 1)1 


A 

3F 

A 

0Q 


AQ** + (1 - a) 


A 

dY 

A 

0Q 


A ** 

AQ/ + 1 


7+1 


(1+0 2 )(A^ 


■ [(j5 - 1 + fj ) n sj - , aq” , - x + 2£ +4 + A* A 9T + ifr + j 5 + .) V+ 1 A 9 T+ i] = 


AQ 


(7.5b) 


Sweep 3 (C direction) 
0i At 


AQI + 


(1 + 0 2 )A£ 


A 

0G 

A 

dQ 


AQ^_ 1 + (2a — 1)( ) AQ2 + (l-a) 


k - 1 


A 

5Q 


A 

0G 

A 

5Q 


A Ql + x 


k + 1 


0 j At 

(1 + 0 2 )(A£) 2 


• [(Tjt— i + fk) n gk- l AQjt,. _ i — ( fk - ] + 2fy, +fk+ i) Sk A Qk + ( fk + fk + l) %k + 1 AQ& + i] 


AQ 


(7.5c) 


The subscripts i, j, and k represent grid point indices in the >j, and C directions. For notational conven- 
ience, terms without an explicitly written ij, or k subscript are understood to be at i,j, or k. In the viscous 
terms on the left hand side, f is the coefficient of djd£ (or <9/ dtj or djd£, depending on the sweep) in the 

dtvJdQ (or aFV,/5Q or dG v Jd Q) Jacobian coefficient matrix. Similarly, g is the term in the parentheses 

following d/0£ (or djdrj or djdQ in the dEyJdQ (or dYyJdQ or dGyJdQ) Jacobian coefficient matrix. 
Equations (7.5a) through (7.5c) represent the three-sweep alternating direction implicit (ADI) algorithm 
used to advance the solution from time level n to n + 1. 
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7.2 MATRIX INVERSION PROCEDURE 
7.2.1 Non-Periodic Boundary Conditions 

The complete set of algebraic equations for the first ADI sweep with non-periodic boundary conditions 
can be written in the following block matrix form. 3 


b'i q 

A 2 B 2 

A 3 

A', 

c 2 

B 3 c 3 




A * 

AQ, 

A * 

aq 2 

A* 

aq 3 

■ 


1 

65* & & ‘ 

1 



• 



• 

= 

• 



■ 

A AT, - 2 B^_2 

C Nj-2 


■ 

aq;._ 2 


■ 

S /V,-2 



A /Vj - 1 

b at,- 1 

C 1V,- 1 

AQaT,-! 


S /V, - 1 



■ C '*l 

A'at, 

B '*, 

aq;, 


S W, 

- 




- 

_ 


- 


These equations result from the application of equation (7.5a) for / = 2 to N\ — 1, with boundary conditions 

A 

added at i » 1 and i=N\. The parameter AQ is the ^-element vector containing the unknown dependent 
variables; A, B, and C are the N eq x N eq coefficient submatrices at /— 1, /, and /+ 1, respectively; and S is 
the ^-element subvector containing the explicit source terms. Also, A', B\ and C' are the coefficient 
submatrices and S' the source teim subvector for the boundary conditions. A variety of boundary condi- 
tions may be used. They are described briefly in Section 6.0, and in greater detail in Volumes 2 and 3. 

Note that the equations at the boundaries may contain coefficients at the boundary point and the two 
adjacent interior points. This occurs, for example, when extrapolation or second order gradient boundary 
conditions are specified. As written, therefore, the coefficient matrix in equation (7.6) is not block 
tridiagonal. However, A{ can be eliminated by multiplying the second row of the matrix by A! Cj 1 and 
subtracting from the first row. CVj can be eliminated in a similar manner. Doing this, we define 

B, = B', — A'j CJ ] A 2 

C, = C', -A', Cj ^ (7.7) 

S^S'j - A', CJ % 

and 

A at, = A'#, — CV, A.v, 1 - iB iVl - 1 

B (Vl = B' (V] -C' A . (7.8) 

S,v, = S'jv, ~ Cjvj A N *_ ]S^ _ , 


3 Although this discussion is written for the first ADI sweep, an exactly analogous procedure is followed for the 
second and third sweeps. 
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The set of algebraic equations solved during the first ADI sweep can now be written as 



A* 


" 

B, C, 

AQ, 


Si 


A« 



A 2 ®2 Q 

aq 2 


S2 

A3 B 3 C 3 

A * 

aq 3 


S3 

■ 

■ 


■ 

■ 

■ 

= 

■ 

■ 

■ 


■ 


A* 



Ay, - 2 - 2 Cy, - 2 

AQat, -2 


Sy, -2 

A at, - 1 ®y, - 1 c y, - 1 

AQy, - 1 


Sy,-i 

A y, By, 

A » 

AQat, 


Sy, 

- 

. 


- 


Since the coefficient matrix is now block tridiagonal, the equations can be solved using the block matrix 
version of the Thomas algorithm (e.g., see Anderson, Tannehill, and Pletcher, 1984). The procedure can 
be summarized as follows: 

1. Define Di = B!. 

2. Compute Ei = Dr ! Ci and AQ5 = Dr ] Si. 

3. For / = 2 to N\, compute 

D; = B; — A ( E, _ 1 

e^d-'q 

AQ; = D“ 1 (S; — A,AQ' _ j ) 

(Actually, E, is only needed for i — 2 to JV, — 1). 

4. Then, set AQ^, = AQ!y, • 

5. Finally, for i = N t — 1 to 1, compute AQ, = AQ! — E,AQ i+ 


In the Proteus code, in step 2 Ei and AQ! are actually obtained by solving DiEi = Cj and D,AQ! — S, 

A 

using LU decomposition of D. A similar procedure is used to compute E, and AQ' in step 3. 

7.2.2 Spatially Periodic Boundary Conditions 

In computational coordinates a spatially periodic boundary 7 condition in the { direction may be repres- 
ented as shown in Figure 7.1. 4 * 


4 As in Section 7.2.1, this discussion is written for the first ADI sweep, but an exactly analogous procedure is followed 

for spatially periodic boundary conditions in the second and third sweeps. 
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Figure 7.1 - Spatially periodic boundary condition. 

The grid points along the / = 1 and i = N\ lines are "similar" in the geometric sense, and have the same 

A A 

flow solution. Therefore, for a spatially periodic boundary condition in the | direction, Qi = Q.V] 

To implement this boundary condition, an additional set of points is added at / = N\ + 1 , setting 

Q Wl + 1 = Q 2 This allows us to use central differencing in the £ direction at i — i'll, computing the coeffi- 
cients in the same way as at the interior points. 
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The resulting set of algebraic equations will consist of N\ — 1 equations (for / = 2 to N\), with N\ + 1 
unknowns. The block coefficient matrix thus has N\ — 1 rows and N\ + 1 columns, as follows: 



A* 



AQ, 


- 


A, 



A2 ^2 

aq 2 


$2 

A3 B 3 c 3 

A * 

aq 3 


S 3 

A4 B4 C4 

A * 

AQ 4 


s 4 

m 

m 


■ 

■ 

■ 

■ 

m 


■ 


A* 



A ty - 2 B ty - 2 ^ty - 2 

A Qty - 2 


Sty- 2 

A ty - 1 B ty - I ^-JV, - 1 

aq;,_. 


Sty - 1 

A ty b at, Cat, 

AQty 


Stf \ 

_ 

aq ;, + 1 




. 



These equations result from the application of equation (7.5a) for /= 2 to N\. As in the previous section, 

parameter AQ* is the jV, 9 -element vector containing the unknown dependent variables; A, B, and C are the 
N eq x N eq coefficient submatrices at /- 1, /, and / + 1, respectively; and S is the jV f9 -element subvector con- 
taining the explicit source terms. 

Since Qj = and Q 2 = Qni + u equation (7.10) can be rewritten with N\ — 1 unknowns as: 


r 

A * 



B2 c 2 a 2 

aQ 2 


^2 

A 3 B 3 C3 

A * 

AQ3 


s 3 

A4 B4 C4 
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A * 

aq 4 

■ 


s 4 

■ 

• 

■ 

■ 

m 


■ 

■ 


A * 



A ty - 2 B ty - 2 C ty - 2 

A Qty-2 


Sty -2 

A ty - 1 B ty - 1 ^AT, - 1 

A * 

AQ v, - 1 

A * 


Sty _ 1 

C.v, a at, b at, 

AQty 


Sty 

_ 

- 


- 


An efficient algorithm to solve this system can be derived that is similar to the Thomas algorithm for 
block tridiagonal systems. The procedure can be summarized as follows: 

1 . Define D 2 = B 2 and F 2 = C*,. 

2. Compute E 2 = D 2 1 C 2 , G 2 = Dj , Aa > and AQ 2 = D 2 ] S 2 . 
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3. For / a* 3 to iVj — I, compute 
Dj = B, — A ; E ; - _ , 

E ; = D-'C f 
Fj = — F / - , Ej _ i 
G ( = -D, ’ Afi; _ , 

AQ; = D~ ! (S ( - — A,AQ' _ , ) 

4. Compute 

G jV, - l ~ D ;Vl ] - i(Cy, _ i — _ jG^ _ 2 ) 

F,v, - i = A, v , - F iV[ _ 2 E, n - _ 2 

D.v, = Bjv, - ^ FA 
i = 2 

AQV, =D^^- £ F,aq' j 

5. Then, set AQ^ = AQVi * 

A A A 

6. Compute AQ Al _ i = AQV t - 1 - G^, _ l AQ A r 1 . 

7. Finally, for / = N\ — 2 to 2, compute AQ, = AQ,- — E,AQ, + j — G,AQ, V , 

* 

In the Proteus code, in step 2 Ej, G 2 , and AQ 2 are actually obtained by solving D 2 E 2 = C 2 , 

A 

\j 2 G 2 = A 2 , and D 2 AQ 2 = S 2 using LU decomposition of D. A similar procedure is used to compute E„ 

A * 

G„ and AQJ in step 3, and _ 1 and AQVi in step 4. 

7.3 UPDATING BOUNDARY VALUES 


With the ADI algorithm described in Section 7.1, if gradient or extrapolation boundary conditions are 
used for the first or second sweep, the boundary values from the first two sweeps must be updated sifter the 
third sweep. This point is easiest to illustrate by looking at a figure. 

In Figure 7.2, a 4 x 4 x 4 grid is shown in computational space for a three-dimensional problem. This 
example assumes that no spatially periodic boundary conditions are being used. The circles represent grid 

A 

points at which the intermediate values Q are computed during the first ADI sweep. These include the 
boundary points at £ = 0 and £ = 1. The squares represent grid points at which the intermediate values 

Q** are computed during the second ADI sweep, including the boundary points at rj — 0 and >7 = 1 . The 

A 

triangles represent grid points at which the final values Q" + 1 are computed during the third ADI sweep, 
including the boundary points at £ ~ 0 and C = 1- If gradient or extrapolation boundary conditions are used 
during the first and/or second sweep, so that the boundary values depend on the interior values, then the 
intermediate values at the £ and/or r\ boundaries must be updated after the third sweep to be consistent 
with the final values at the interior points. 

To do this, after the last sweep the difference equations are rewritten and solved at the <£ and y\ bound- 
aries. At the £ = 0 boundary, 
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B\ *AQ, + C\ n AQ2 + A\ "AQ" = S', 71 (7. 12) 

The subscripts refer to the value of z, the index in the £ direction. This equation is applied for j — 2 to 
A /2 — 1 in the rj direction, and for k = 2 to N$ — 1 in the £ direction. For notational convenience, however, 
the subscripts j and k have been omitted. 


All the terms in equation (7.12) are known except AQ?. Solving, 


AQ” = (B'i V ’(S', " - C\ "AQ” - A', n AQ”) 


At the £ = 1 boundary, 


Cat, _ 2 + AV, "AQ^ _ , + B' A) AQ^ = S' A - " 

= (B^ n ) \S' Nx n - CV t "AQ^ _ 2 — A'^ "AQ^ _ ,) 

An analogous procedure is followed to update values at the rf = 0 and rj = 1 boundaries. 


(7.13) 

(7.14) 

(7.15) 


Finally, note from Figure 7.2 that new values are not computed at the comers or edges of the compu- 
tational domain during the solution algorithm. To make the edge values consistent with the rest of the flow 
field, in Proteus they are defined using the computed values at adjacent points. For example, at each £ lo- 
cation along a £->7 edge (i.e., one of the four lines of intersection between the £ and rj boundary planes), the 
density p and total energy E r are arbitrarily defined by linearly extrapolating from the two adjacent points 
in the £ and r\ coordinate directions, and averaging the two results. The edge values of the velocities are 
updated by doing the same type of extrapolation. Instead of averaging, however, the extrapolated velocity 
whose absolute value is lower is used. This was done to maintain no-slip conditions at duct inlets and exits. 
The values at comers, where all three boundary planes intersect, are determined in an exactly analogous 
manner. 
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Updating boundary values from the first two sweeps is complicated somewhat when spatially periodic 
boundary conditions are used. Details are presented in the description of subroutine BVUP in Volume 3. 
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8.0 ARTIFICIAL VISCOSITY 


With the numerical algorithm of Section 7.0, high frequency nonlinear instabilities can appear as the 
solution develops. For example, in high Reynolds number flows oscillations can result from the odd-even 
decoupling inherent in the use of second order central differencing for the inviscid terms. In addition, 
physical phenomena such as shock waves can cause instabilities when they are captured by the finite dif- 
ference algorithm. Artificial viscosity, or smoothing, is normally added to the solution algoritlun to suppress 
these high frequency instabilities. Two artificial viscosity models are currently available in the Proteus 
computer code - a constant coefficient model used by Steger (1978), and the nonlinear coefficient model 
of Jameson, Schmidt, and Turkel (1981). The implementation of these models in generalized 
nonortho gonal coordinates is described by Pulliam (1986b). 

8.1 CONSTANT COEFFICIENT ARTIFICIAL VISCOSITY 

The constant coefficient model uses a combination of explicit and implicit artificial viscosity. The 
standard explicit smoothing uses fourth order differences, and damps the high frequency nonlinear insta- 
bilities. Second order explicit smoothing, while not used by Steger or Pulliam, is also available in Proteus . 
It provides more smoothing than the fourth order smoothing but introduces a larger error, and is therefore 
not used as often. The implicit smoothing is second order and is intended to extend the linear stability 
bound of the fourth order explicit smoothing. 

The explicit artificial viscosity is implemented in the numerical algorithm by adding the following terms 
to the right hand side of equation (7.5a) (i.e., the source term for the first ADI sweep.) 

(V^Q + V^Q + V { A { Q) - [(V ? A f ) 2 Q + (V^Q + (V ; A ? ) 2 Q] (8.1) 

where and tf are the second- and fourth-order explicit artificial viscosity coefficients. The symbols V 
and A are backward and forward fust difference operators. Thus, 

— Q i ~ Q/ - 1 
i — Qi + 1 — Q i 

VjQ, = Q/ + 1 ■— 2Q, + Q/-] 

(V ? A { ) 2 Q f = Q i + 2 - 4Q; + , + 6Q, -4Q,_, + Q;_ 2 

Equivalent formulas are used for differences in the >/ and C directions. 

A few details should be noted at this point. First, the sign in front of the artificial viscosity term being 
added to equation (7.5a) depends on the sign of the "i" term in the difference formula. For damping, that 
term must be negative when added to the right hand side of the equations (i.e., explicit artificial viscosity), 
and positive when added to the left hand side (i.e., implicit artificial viscosity.) See Anderson, Tannehill, 
and Pletcher (1984) for details. Second, the terms being added are differences only, and not finite difference 
approximations to derivatives. They are therefore not divided by A£, etc. Third, the variables being dif- 
ferenced are Q, not Q. As noted by Pulliam (1986b), scaling the artificial viscosity terms by 1/7 makes them 
consistent with the form of the remaining terms in the equations. Fourth, the terms are also scaled by At. 
This makes the steady state solution independent of the time step size (Pulliam, 1986b). And finally, note 
that the fourth-order difference formula cannot be used at grid points adjacent to boundaries. At these 
points, therefore, the appropriate fourth-order term in expression (8.1) is replaced by a second order term. 
Thus, for points adjacent to the £ = 0 and £ = 1 boundaries, — £? ) At[(V { A { ) 2 Q]/7 is replaced by 
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A similar expression is used at points adjacent to the rj and C boundaries. 

The implicit artificial viscosity is implemented by adding the following terms to the left hand side of the 
equations specified. 


£ /At a, 

-L-[V^(JAQ )] 
£ /At a** 

)] 

-^[V C A C (/AQ")] 


to equation (7.5a) 
to equation (7.5b) 
to equation (7.5c) 


(8.3) 


Note that the addition of the artificial viscosity terms, in effect, changes the original governing partial 
differential equations. At steady state, the difference equations with the artificial viscosity terms added ac- 
tually correspond to the following differential equations. 5 


dE dF dG 

d£ dr, + 5C 


A A 

dE v d¥ 

■ + 


dr, 


A 

y dG 

+ 




f (2) 
V , l E 

J 


2 d 2 (JQ) 2 d 2 (JQ) 2 

(A £) 2 — + (A>j ) 2 — + (AQ 2 


d? 


drf 


J 


4 d 4 (7Q) 4 d 4 (JQ) „ 4 d 4 (JQ) 

(AC ) 4 — ^ + (A ^) 4 — + (AC ) 4 


ar 


dn 


dC 


9 A 

d 2 (7Q) 

dC 2 


A 

The implicit terms do not appear, since they difference AQ, and in the steady form of the equations 

AQ = 0. The artificial viscosity terms do not represent anything physical. The coefficients should therefore 
be as small as possible, but still large enough to damp any instabilities. Although optimum values will vary 
f rom problem to problem, recommended levels are = C>(1) and £/ = 2sf> (Pulfi^m, 1986b). The recom- 
mended level for when used, is £?> = 0(1). 

8.2 NONLINEAR COEFFICIENT ARTIFICIAL VISCOSITY 


The nonlinear coefficient artificial viscosity model is strictly explicit. Using the model as described by 
Pulliam (1986b), but in the current notation, the following terms are added to the right hand side of 
equation (7.5a). 


+ V w 


+ V, 


( “7 ) + ( "j ) ( c ? : >A {Q *" 4 4>A £ v r A {Q)« | 

(t) + (4)l(*?V-4 4 VAO)y; 

+ 1 v / j mm 

(t) +(t)Vv-4 4 Vaq).| 

^ ' k + 1 ^ h 


(8.4) 


5 These equations represent the use of the constant coefficient artificial viscosity model presented in this section. The 
nonlinear coefficient model to be presented in Section 8.2 is more complicated, but the same principle applies. 
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The difference operation A^V^A^Q is given by 

AjVjAjQi* = Q/ + 2 “ i + 1 ^ Q/ — 1 

In the expression (8.4), \j/ is defined as 

^ = tx + ty + 

where \J/ X , i// y and \j/ 2 are spectral radii defined by 6 

\u\ +ay[ZF+tf^zF 


(8.5) 


'l'x = 


A{ 


^v = - 




\V\+aJ 


2 2 2 
r\ x +V y +1 Z 


( 8 . 6 ) 


A»j 


| W\ + a v 'C * 2 + C y 2 + C 2 2 
AC 


Here 6/, F, and IF are the contravaiiant velocities without metric normalization, defined by 

U = i t + £ x u + C y v + 

V = It + njc w + V + (8.7) 

TF=C, + C^ + Cj,v + C 2 w 

and a = , JyRT , the speed of sound. 

The parameters t Q) and are the second- and fourth-order artificial viscosity coefficients. Instead of 
being specified directly by the user, as they are in the constant coefficient model, in the nonlinear coefficient 
model they are a function of the pressure field. For the coefficients of the £ direction differences, 


where 


(e| 2) ), = k 2 At max(cr i - + , , a v a t _ J 
(4 4) )/ = max[0, jc 4 At - (cf } )J 

_ Pi + 1 ~ jft + Pi - 1 
P i+ i+2Pi+Pi-i 


(8.8a) 

(8.8b) 

(8.9) 


Similar formulas are used for the coefficients of the rj and C direction differences. 

The parameter c is a pressure gradient scaling parameter that increases the amount of second order 
smoothing relative to fourth order smoothing near shock waves. The logic used to compute s w switches 
off the fourth order smoothing when the second order smoothing term is large. 

The parameters tc 2 and k 4 are user-specified constants. Like the coefficients in the constant coefficient 
model, the optimum values will be problem-dependent, and are best chosen through experience. Cases have 
been run with values of k- ranging from from 0.01 for flows without shocks to 0.1 for flows with shocks, 
and k 4 ranging from 0.0002 for flows computed with spatially constant second-order time differencing to 


6 It should be noted that the grid increments A£, A r\, and AC in these definitions do not appear in the corresponding 
formulas presented by Pulliam (1986b). This is because the grids used by Pulliam are constructed such that 
A£ = Arj = A{ = 1, while in Proteus AC = l/(Ni - 1), Atj = 1 l(N 2 - 1), and AC = l/(A r j - 1). The definitions used 
here for \j/ x , v y , and i j/ 2 result in an artificial viscosity level equivalent to that described by Pulliam. 
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0.005 for flows computed with spatially varying first-order time differencing. Pulliam (1986b) gives 
k 2 ~ 0.25 and k 4 = 0.01 as typical values for an Euler analysis. 

Like the constant coefficient artificial viscosity model, the nonlinear coefficient model requires special 
formulas near boundaries. To apply (8.4) at i = 2, cf> is needed at / = 1. It is defined as 

(e^)j = k^At max(<72, CT i) 

With the above definition, applying (8.4) at / = 2 and / = M — 1 requires a at i = 1 and i = AV They are 
defined as 

— Pa + 4Pi - 5Pi + 2 P\ 

° X ~ P4 + 4P3 + 5 P2 + 2 P\ 

- Pn x - 3 + 4P,V, - 2 ~ 5 />V, - 1 + 2 PN X 
Pn x - 3 + 4 Pn x - 2 + 5 P/f x - 1 + 2 Pn x 

And, finally, applying (8.4) at / = 2 and i = N\ - 1 requires A^A^Q at i = 1 and i= N\— 1. There are 
numerous formulas that could be used. The ones currently in the Proteus code are 

AjVjAjQj = — Q 5 + 5Q 4 — 9Q 3 + 7Q 2 — 2Qj 
A^V^Q^i _ j = Q N[ _ 4 — _ 3 + 9Q^ _ 2 “ - 1 + 
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9.0 TURBULENCE MODEL 


As noted briefly in Section 2.0, for turbulent flow the Reynolds stress and turbulent heat flux terms are 
modeled using the Boussinesq approach. An effective viscosity is thus defined as fi — pa + pit, where & is 
the laminar, or molecular, viscosity coefficient, and \i t is the turbulent viscosity coefficient. Similarly, an 
effective second coefficient of viscosity is defmed as X = Xt + X : , and an effective thermal conductivity coef- 
ficient is defined as k = k t + K 

The turbulent coefficients must be computed using a turbulence model appropriate for the flow being 
computed. In Proteus, turbulence is modeled using either a generalized version of the Baldwin and Lomax 
(1978) algebraic eddy viscosity model, or the Chien (1982) low Reynolds number k-z model. 

9.1 BALD WTN-LOMAX MODEL 

For wall-bounded flows, (i.e., boundary layers), the Baldwin-Lomax turbulence model is a two-layer 
model, with 


[(dinner for y n <y b 

Hr=< (9- 1 ) 

(_ (Mf) outer for y n >y b 

where y n is the normal distance from the wall, and y b is the smallest value of y n at which the values of \x t from 
the inner and outer region formulas are equal. For free turbulent flows (i.e., mixing layers, jets, and wakes), 
= (fj, t ) ouUr . In the inner region, in addition to the Baldwin-Lomax model, an alternate expression first 
presented by Spalding (1961), and later by Kleinstein (1967), is also available. 

9.1.1 Outer Region 

The outer region turbulent viscosity is computed from 

(M/) outer = KCcppFfciefrF , vdkeRCr (9*2) 

where K is the Clauser constant, taken as 0.0168, C cp is a constant taken as 1.6, and p is the static density. 
The parameter F wake is computed from 


Fwake 


| ymaxFmax 
n t/ 2 ymax 

wk * diff p 

r mnY 


for wall-bounded flows 
for free turbulent flows 


where C wk is a constant taken as 0.25, and 


v w= v 


V 


(9.3) 


where V is the total velocity vector. 

The parameter F max in equation (9.3) is the maximum value of 
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F(jn) = 


y* |«| 


for wall-bounded flows 
for free turbulent flows 


(9.4) 


and y mex is the value of y„ corresponding to 

In a simple boundary layer analysis, with only one solid surface, the procedure for computing F(y n ) and 
F max is relatively straightforward. In a general Navier-Stokes analysis, in which any part of any boundary 
may be a solid surface, the problem is more complicated. 

In Proteus , each grid point is labeled as either a wall-bounded point or a wake point. Wall-bounded 
points are those for which at least one of the three grid lines through the point intersects a solid wall. All 
other grid points are wake points. 

For each wall-bounded point P wdh the grid line intersecting the nearest wall is determined. F(y n ) is then 
computed along that line, with y n equal to the distance to the wall nearest the point P wct! . F max is defined 
as the maximum value of F along the line, and y max is the value of y n corresponding to F mgx . It has been 
found that for wall-bounded flows the function F(y n ) can have multiple peaks. In 3-D Proteus , the peak 
nearest the wall is used. 

■ For each wake point P wa k ei the grid line intersecting the nearest boundary is determined. Next, the values 
of I V\ and \V\ on that line are found. Two values of F(y n ) are then computed - one with y n equal 
to the distance from the point P wakt to the location of | V\ . , and one with y n equal to the distance to the 

location of | V\ . Two values of F max and y max are determined, for the two F(y n ) arrays. As in the wall- 
bounded case, yZll is the value of y n corresponding to F max . The smaller and the corresponding F max are 
the values finally used for computing F wak e . 


In equation (9.4), Q is the magnitude of the total vorticity, defined as 



/ dw dv \ ( du dw \ 2 / 

\ dy dz ) \ dz dx ) \ 


dv 

dx 



(9.5) 


The parameter A + is the Van Driest damping constant, taken as 26.0. The coordinate y+ is defined as 


-§- Pw^yn „ yJ r wP w^ e r 

y — — **= — s — y* 


where u, = Jr„lp„Re r is the friction velocity, r is the shear stress, and the subscript w indicates a wall value. 
In Proteus, t „ is set equal to | Q ] . 


The function F K m in equation (9.2) is the Klebanoff intermittency factor. For wake points, F Kleb = 1. 
For wall-bounded points, 


r Kleb {^K!eb)mm 4 " [1 (^Kteb^mini 


1+5 


( 


^KlebFn 

Fmax 



(9.7) 


This factor accounts for the experimentally observed fact that, as the free stream is approached, the fraction 
of time the flow is turbulent decreases. In equation (9.7), 5 and Ckm are constants taken as 5.5 and 0.3, 
respectively. (C Klcb ) mm is a constant normally equal to 0.0. However, when using the Baldwin-Lomax model 
to generate initial turbulent viscosity values for the Chien k-c model (discussed in Section 9.2), is 

set equal to 0.1. This yields a small positive value for /x, in the free stream, and has been found to minimize 
starting problems with the k-t model. 
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9.1.2 Inner Region 


The inner region turbulent viscosity in the Baldwin-Lomax model is 

(m,w=p/ 2 |oK (9 - 8) 

where / is the mixing length, normally given by 

/-■%(! -e-”' 14 ') ( 9 - 9 ) 

and k is the Von Karman constant, taken as 0.4. 

A modified form of equation (9.9), proposed by Launder and Priddin (1973), may also be used. This 
formula is most useful for flows with steep negative gradients of shear stress normal to the wall, such as 
accelerated flows or flows with suction. Their modified formula for / is 

l= K y n (l-e- y ^ y!A+ ) (9.10) 


where 


+ 


T 



and n is a constant taken as 1.7. 

The inner region turbulent viscosity may also be computed using an alternate expression first presented 
by Spalding (1961), and later by Kleinstein (1967). In this model, 

(fiiimtr - -l-Ki/ + --±- (,» + ) 2 ] (9 11) 

where 


T \! r wl Pw^ € r 


Again, in Proteus, t* is set equal to | £2 1 . 

9.1.3 Turbulent Values of A and k 

The turbulent second coefficient of viscosity is simply defined as 

(9 12 ) 


The turbulent thermal conductivity coefficient is defined using Reynolds analogy as 


‘ Pr t r 


(9.13) 


where c p is the specific heat at constant pressure, and Pr t is the turbulent Prandtl number. In Proteus , the 
turbulent Prandtl number may be treated as constant, or as a variable using the following formula (Wassel 
and Catton, 1973): 
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Pn= 


c, 


I — exp[ — 


Pr3 


CprA \ 
I 1 tit 1 ! ) 


C Pn Pr, 


1 — exp( — 


C, 


PrinJm ) 


(9.14) 


Here C Pr \> C Fr i, C P and C P * are constants taken as 0.21, 5.25, 0.20, and 5.0, respectively, and Pn= c^ki 
is the laminar Prandtl number. 

9.2 CHIEN krt TURBULENCE MODEL 

9.2.1 k-z Equations 

The low Reynolds number k-z formulation of K. Y. Chien (1982) was chosen because of its reasonable 
approximation of the near wall region and because of its numerical stability. Here k and z axe the turbulent 
kinetic energy and the turbulent dissipation rate, respectively. 7 In addition, the Chien k-e turbulence model 
was frequently used in past Navier-Stokes computations with good results (Nichols, 1990, 1991; Patel, Rodi, 
and Scheuerer, 1985; Sahu, 1984.) The set of k-z equations are lagged in time and solved separately from 
the Navier-Stokes equations to allow for code modularity in turbulence modeling. In Cartesian coordinates, 
the three-dimensional equations for the Chien k-z model can be written using vector notation as 


aw , aF , aG , aH „ . ^ 

—7; — b “5 ‘ 5 t — “ — - of 1 

dt dx dy dz 


(9.15) 


where 


F = 


G = 


H = 


S = 


w = 

{*} 

puk — 

Re r ^ 


1 

puz — 

Re r 

pvk — 

Re r ^ 


1 

pvz — 

Re r ^ ' 

pwk — 

Re r ** 


1 

pwz — 

Re r ^ 

p k 

— Re r pz 


dk 

dx 

dz 

dx 

dk ' 

dy 

dz 

dy 

dk 

dz 

dz 

dz 


C } P k f-Re r C 2 p i- 


(9.16a) 


(9.16b) 


(9.16c) 


(9. 16d) 


(9.16e) 


It should be noted that in the Chien model, t is actually the isotropic portion of the turbulent dissipation rate. 
Throughout this manual, however, it is referred to as simply the turbulent dissipation rate. 
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2 M 


T = 


Re r 

2 Mg 
Re r 


y n 

-y + ! 2 


2 

J'ii 


(9. 16f) 


and 


. Mr 

(9.17a) 

, Mr 

m £ = m + ^- 

(9.17b) 

C, = 1.35 

(9.17c) 

C 2 = Q,( l-|e"^ /36 ) 

(9.17d) 

°k =i-° 

(9.17e) 

^=1-3 

(9. 1 7f) 

oo 

II 

(9. 17g) 

R P k \ 

K t ~ pi 

(9. 17h) 

^-Re, P '- 

(9.18a) 


/>,=2 


/ 5m \ 2 . / dv \ 2 ( dw \ 2 2 / du , dv , dw \ 

\ dx ) \ dy ) \ dz ) 3 \ dx dy dz ) 
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dv 
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dw 

dw 

— + 

— + 
dz + 

dx 
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— du 
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dw 



*2 

dx 

! 

dy 
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The turbulent viscosity is given by 


H t =C u p±- 


c ,=9>.(‘ 




C u =0.09 

rr 


C 3 = 0.0115 


(9.18b) 

(9.18c) 

(9.19a) 

(9.19b) 

(9.19c) 

(9.19d) 


Note that the vectors VV, F, G, H, and S are used in most standard k-z. formulations (with different 
constants), and the vector T is unique to the low Reynolds number formulation of Chien. The parameter 
y n is the minimum distance to the nearest solid surface, and y + is computed from y„. The production of 
turbulent kinetic energy P k includes the full Boussinesq approximation for compressible flows. All of the 
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above equations have been nondimensionalized using appropriate normalizing conditions. 
Nondimensionalization of mean flow properties is discussed in Section 2.1. The turbulent kinetic energy 
k and the turbulent dissipation rate z have been nondimensionalized by u } and p r itflti r , respectively. 


Following the procedure of Section 2.3, the following generalized grid transformation is used to trans- 
form the k-z equations from physical (x, y, z, t) coordinates to computational (<*, f, t) coordinates. 


£ = S(x,y f z) 

*i = v(x>y>z) 

C = C(x,y,z) 

T= t 

Applying the generalized grid transformation to equation (9.15) yields 


(9.20) 


VV + + F ?i x + F Z x + G*£ y + G n rj v 


i r* r . u 


(9.21) 


Although the above equations can not be put into exact strong conservation law form, the procedure 
used to do so for the mean flow equations, described in Section 2.4, is nonetheless applied to equation 
(9.21). The result is 


where 


0 W. 3 F.aG.dH £ . £ 

~dT + ~dT + ~d^ + ~dt'~ i>+ 1 


w= — 

J P* J 


A A 

F = F, 


A A 

F n - F. 


A 1 

Fc= y 


C~ r D~ r M 

i x puk + Zypvk + £ z pwk 
Z x pu£. + £ y pvz + Z z pw£ 


A 

Fn = 


J 1_ 

J Re r 


t*M+Zy+& k t 


(9.22) 

(9.23a) 

(9.23b) 

(9.23c) 

(9.23d) 


A 



J 1_ 

J Re r 


Pi kiZx’lx + Zyty + t&JKi - Pk(Zj£x + Zy£y + 
Peitxlx + Zyly + ~ PlUx + l£y + 


(9.23e) 


G c - 

A 

G Z) = 


AAA A 

G — G r — Gn — G 


M 


r\ x puk + r]ypvk + rj 2 pwk 
y] x pus + Y{ypvz + 


1 1 


J Re r 


Pkinl + >l 2 y + 1 2 z) k v 

, 2 , 2 , 2. 

1 1 I v? - 4 — v? 4 - vi k 


(9.23f) 

(9.23g) 

(9.23h) 


'X 11 Pki^xtfx "F %y*ly "F Pki*! x£x "F "F 

^ J Re r Pzi^xtfx "F ^y*ly "F — Pcfax^x "F *ly£y "F 


(9.23i) 


H = H C -H i) -H w 


(9.23j) 
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A 1 
Hc "T 


C xPuk + C ypvk + C z pwk 
£ x pue + C y pV£ + C^pwe 


1 1 


Hp J Re r 


^l + fy + tik 


A 1 1 
H M = - 


J Re r 


Pk^x^x £y£y ^z^z)^l Pki^x^x ^z^z^r, 

P £ (c .xCx + CyCy "b ^ z^r) f^ti^x^x "b *ly^y "b r iJr.z) t r, 


S = - 


P k - Re r pz 

C,P k -j-Re r C lP \ 


(9.23k) 

(9.231) 

(9.23m) 

(9.23n) 


T = t 


2 P 


Re r -2 


2 P e 


y n 

-y + R 


Re r 


2 

y n 


(9.23o) 


Note that in equation (9.23n), the term P,. involves derivatives with respect to the Cartesian coordinate 
directions (see equations (9. 18a-c.) These are evaluated using the chain rule. 

9.2.2 Linearization of the k-z Equations 

Solving equation (9.22) for 3W/3 t and substituting the result into the time differencing scheme of Beam 
and Wanning (1978), given by equation (3.1), for d(AW")/dT and dW"/d t yields 


AW": 


1 + 02 \ 

dl dl 

■ 

<T> 

+ 

1 

drj dr} 

drj 

dl 

1 A7 

( 3 F c . df D . 

dc c 

dG D dG M 

dHr 3 H C 


i+a 2 

a? aj 

dl dr, 

dr} drj 

dl ' 81 ' 

dl 

+ * 

+ i +e 2 

aw* -1 + 0 ^( 0 ,-- 

±-e 2 )(Ax)\ 

(At) 3 ] 




dl 


dl 


- + AS + AT 


Equation (9.24) is then linearized using the procedure described in Section 4.0. Let 


B = ^-, C = .. 

aw aw aw 


d* D ^ d= ac 0 , r _ 3 H c f = y.JT. 


aw 


aw 


aw 


aw 


aw 


be the Jacobian coefficient matrices, where 

A = 


l x U + {j,V + Z z W 

0 


0 

ixU + ZyV + £ 


J 


( 9 . 24 ) 


( 9 . 25 ) 


(9.26) 


B = 




JRe r <*<<d + <l + il)( p ) 


(9.27) 
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c= 


rj x u + ri y v + rj z W 0 

0 Y] x u + rjyV + rj z w 


(9.28) 


D — 


JRe r Ski'll + t i + r l 2 j{ p ) 


JRe r + r, y + P ) 


(9.29) 


£ = 


+ C Z * 

0 


0 

+ C v v + C Z W 


] 


(9.30) 


F = 


JRe r *k(£ + Z 2 y + &( p ) 


JRe r + + P ) 


(9.31) 


M = 


or 1 A A 

Z< -> £ M f 


/>* s 2 

Q^-^+Pe r C 2 ^- 


k 2 ?k 

— C — — — Re 

S* 2 Mf ' 
£ 

-2Re r C 2 f 


(9.32) 




2m 


pyl Re r 




py 2 n Re r 


(9.33) 


The linearized form of equation (9.24) can now be written as 

, u > . r 3 < gAvV > , , AgAl _ _ v ( a^T = 

+ i + e 2 I a^ dr/ ac ac J 

0 ,At I" a(AF M ) a(AG w ) a(AH w ) “I 

i+a 2 |_ a? + dr, a? J 

/ aF c dFn aG c aG 0 aG M aH c aH^ aA^ * »\ 

^ d* + dl d% drj drj drj d£ d£ d£ J 


At 


M , c ■ T 1 i ^ AW^ 71 - 


i +e 7 


+ o[(a,-^-a 2 ) ( At )2 , ( A' )3 ] 


( 9 . 34 ) 
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9.2.3 LU Factorization Algorithm for the A-e Equations 

The LU factorization scheme used to solve the k-s equations is essentially the same as that described 
by Hoffmann (1989). Letting RHS(9.34) represent the right hand side of equation (9.34), we can write 


1 + 




RHS(9.34) (9.35) 


where I represents the identity matrix. 

The Jacobian matrices A , C, and E can each be split into two submatrices, such that each submatrix is 
associated with the positive and negative eigenvalues of the corresponding Jacobian matrix. Equation (9.35) 
can thus be rewritten as 

{ ,+ T^[w < ' , * + - , ' )+ i (C * +c ' ,+ i <£ * +£ " > -'if-f'''t- <M+w ]}" 4 ' v * RHS|9 - 3 ‘ ) <,!6) 

Using first-order upwind differencing for the Jacobian matrices A, C, and E , and central differencing for 
B , D , F , and on the right hand side, equation (9.36) becomes 

+ \ l +e 2 + + 6 < A ~ + 6 n C+ + 6 n C ~ + 6 i E+ + S t E ] 

- Trk + (6 < Br + {6 * D)+ + (6 ' Dy + (<5 « f)+ + (6 < n ~ +{M+ ^ = i+k + ’ 1 

+ j At ^ ( - <5 { F C + + “^F^ — 4,G C + <5,G C + 6 q G m — <5 { H C + 6 ( H D + + S + t) (9-37) 

+ YTe; AW" - 1 + o[(^, - \ - e 2 )(Ar) 2 , <aV) 3 J 

Note that the central differencing operators for B } D, and F have been split into forward and backward 
differencing parts. Neglecting the temporal truncation and splitting errors, equation (9.37) can be approxi- 
mately factored as 

-v n 

j 1 + ^ 6 < A + + ^ c+ + 6 < E * ~ {6<Bf ~ {6 ^ r ~ <vnj • 

jl + [ [6+A- + 6+C~ + <5j £ - - (6 ( B) + - (<5,D) + - (6 ( F) + — (M + /V)]| AW" = RHS(9.34) [<5 { (AF W ) + <5, (AG*) + <5 ( (AH*)] n 


At 

1 + 0. 


- ( — <5^F C + <5^F C + <5 { F* — <5,G C + 6 n G D + <5,G* — <5 { H C + <5 f H 0 + <5 f H* + S + T) 


(9.38) 


° 2 AW"’ 1 


1 +6 


This equation can then be split into the following two-sweep sequence. 
Sweep 1 (upward) 


jl + [ SJA + + 6~C+ + 6^E + - {6,8)- - {5 V D)~ - (^/)"]| AW* = RHS(9.34) (9.39a) 


Sweep 2 (downward) 


|l + [si A- + S+C- + S* ( E- - ( S t B ) + - (i,D)+ - (V0 + ]| AW" = AW‘ 


(9.39b) 
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9.2.4 LU Sweeping Procedure for the k~c Equations 


Non-Periodic Boundary Conditions 

In the solution algorithm, the upward sweep is done first, then the downward sweep. When applied at 
an interior point in the computational domain, equations (9.39a) and (9,39b) each have just one unknown. 
These equations are therefore solved point-by-point. 

The upward sweep starts at the lower left front comer of the computational domain, point (2, 2, 2), and 
marches along planes of constant i+j+k to the upper right back comer, point (i\\ — 1, N 2 — 1, iV 3 — 1). 

Equation (9.39a) is solved for the intermediate unknown AW* at each point (/,/, k) using known informa- 
tion at points (/— l t j,k) f ( ij — 1 y k), and ( ij\k — 1). This is possible because the left hand side of the 
equation contains only backward differencing operators. 

The downward sweep is in the opposite direction, from point (N\ — 1, N 2 — 1, N 2 — 1) to point (2, 2, 

2), again along planes of constant / +j + k. Equation (9.39b) is solved for the final unknown AW* at each 
point (ij, k) using known information at points (f + 1 J, k ), (ij + 1, k ), and (ij, k -I- 1). This is possible 
because the left hand side of the equation contains only forward differencing operators. 

Spatially Periodic Boundary Conditions 

A spatially periodic boundary condition in the £ direction may be represented as shown in Figure 7.1. 

A A 

Following Section 7.2.2, an additional set of grid points is added at /= N\ + 1, setting W^ + i = W 2 . This 
allows us to use central differencing in the £ direction at / = N\. The upward sweep therefore goes from 
point (2, 2, 2) to point (N u N 2 — 1, A 3 — 1), and the downward sweep goes from point (N u N 2 — 1, 
iV 3 — 1) to point (2, 2, 2). An analogous procedure is used for the periodic boundary conditions in the r\ 
and/or £ directions. 

9.2.5 Updating Boundary Values for k-z Equations 

For easy modification and easy accommodation of complicated boundary conditions for k and e, non- 
periodic boundary conditions are treated explicitly in the solver. After the k and t values at the interior 
points are advanced in time, the values at the boundaries are simply computed from the new interior values 
using the specified boundary conditions. 

Spatially periodic boundary conditions in any sweep direction are treated implicitly, as described in the 
previous section. For a periodic boundary condition in the £ direction, the k and t values at / = 1 are easily 

A A 

updated by setting W t = W^. An analogous procedure is used for periodic boundary conditions in the tj 
and £ directions. 

9.2.6 Turbulent Values of A and k 

The turbulent second coefficient of viscosity X t and the turbulent thermal conductivity coefficient k t are 
defined as described previously in Section 9.1.3. 
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APPENDIX A - EXPANSION OF VISCOUS TERMS 


In Section 4.2, the viscous terms in the governing equations are linearized. To do this, the elements of 

E k , FV, and Gy, given in equations (2.17e) through (2.17g) must first be rewritten in terms of the dependent 
variables, and with derivatives in the Cartesian directions transformed to derivatives in the computational 


directions using the chain rule. The non-cross derivative terms, involving E^, Fk,, and G Vl , are then 

° AAA 


linearized using Taylor series expansion. The cross derivative terms, involving E y 2 , Fp 2 , and G^ 2 , are simply 
lagged one time step. This Appendix presents the fully expanded viscous terms required in the linearization 
procedure. 


The viscous term is given by equation (2.17e), which is repeated here. 




l l_ 

J Re r 


0 

xx 1 * x d" T xy%y T xz^z 


T 
T 

A / 

x x 1 '■yz'^y 

P X Z X + l SyZy+dztz 


xy?x + T yy^y + r yz^z 
xz^x + T yz^y + T zz^z 


(A.l) 


where 


T y-j- = 2flU X + /.{U x + Vy + W z ) 
Tyy ~ 2flV y + ?-{U x + V y + W 2 ) 
t 22 = 2fiw z + X(u x + Vy + w 2 ) 



r xy 

= p(Uy + V x ) 




x xz 

= n{u z + W x ) 




T yz 

= M(V 2 + Wy) 



X 

II 

ux xx + 

VT xy + VVT X2 - 

1 

Pr T 


II 

ux xy + 

VTyy + WTy 2 — 

1 

Pr T 

q y 

Pz = 

« XZ + 

VT yz + WT 22 - 

1 

Pr r 

4z 


<7x=- kT x 

q y = - kT y 
<7z = - kT z 


The chain rule is used to transform derivatives in the Cartesian directions into derivatives in the com- 
putational directions, resulting in 

= (2 y. + + VjH, + + Ktyt + Vn + W + + VSj + C. >{) 
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^yy ( 2 M "f" 2 )(£yVj *ly V ri d” ^(^x^[ ~^~ *lx^rj "b Cx^f) ~f 2(£ 7 Wj “b ^Iz^tj “b C 7 ^(j) 

* 77 = (2 fx + l)(Z 2 W£ + fi 2 w v + t z w Q ) + Klx u > + Vx u r, + Cjc w c) + ; -(Vf + V* + £y v f) 

T xy = M(C y w, + ri^ + Cy« c + ^ ^ + C^) 

J X7 = mCCt^j "b ^Z U rj “b ^Z W £ "b ^X^f "b r lx W t] ~b Of) 

V = m(£ 7 ^ + ^ + C 7 v c + ^Wj + ^ + Cj,w c ) 

Px = ( 2 M + ^)(£x““{ + ’IxUUr, + Cx 1 ^) 

+ MZyUVf + rjyUV^ + CyUVt- + Z 2 UW^ + VzUWr, + Cz^c) 

+ **(£/«{ + *lyVU v + ZyVU(- + £ x Wj + r , x + £ x Wj 

+ m( 5 7 ^ + + C 7 ww c + ^ww ? + ^ + C^mvr) 

+ k ^X T l + Vx^r, + Cx 7 ^) 

Py = (2M + 7 )(^VVj + + C y VV ? ) 

+ 2 (^x v ^ + »?x VM ,, + ^ vw c + + n 7 ^ n + C 7 vw c ) 

+ niZyUUf + rjyUU^ + ZyUUr + ^UVj + rf x UV^ + t x m 0 

+ m(£ 7 *^ + >7 7 wv^ + C 7 wv e + ^HWj + riyWW^ + CyWW^) 

+ k(i y T J + rj y T v + CyTr) 

Pz = (2m + 2)(£ 7 ww f + rtjvw^ + C 7 ww c ) 

+ x wuf + rf x wu n + C x wu q + Zy w s + + Cy^) 

+ n{z^ + y] 2 uu n + C 7 ^ + ■fx^ + Vx^ + Cjc^c) 

+ M(£ 7 W! + *? 7 % + C 7 W f + ?yVWj + IJ^VH^ + CyVW^) 

+ K£ z Tg + >^7"^ + C Z T^) 


The above expressions for the t's and p's are next substituted into equation (A.l). The p derivative 

A A 

terms become elements of Ev p and the tj and £ derivative terms become elements of E v 2 . The resulting five 

A 

elements of Ek[ (excluding the 1 jJRe r coefficient) are 

(Ej/j), = 0 (A. 2a) 

(E k ,) 2 = 2 M ?x^ + 2^x(^f + Zy v t + %z w t) + tfy(£y“S + ^x v j) + + ^x w j) ( A - 2 b) 

(E k ,) 3 = 2 m^v ? + Z£y(Zx“i + £y v i + + + £* v f) + vZz(Zz v i + Z y w i) ( A -2c) 

(Ei/,) 4 = 2 + t£z^x u l + Zy v l + ijVf) + ^x(Zz u l + £x w f) + rfy(Zz v S + {y w f) ( A -2d) 

(E^) 5 = 2/u(£xUU> + ^yWf + c 7 >^’j) + ^x(^x uu l + Zy w l + Zz^t) 

+ ?.£y(Z x vur + ZyWf + ^ 7 wv ? ) + A£ Z (Z X WU^ + ZyWV^ + ^«Wj) 

+ ^x^y vu l + ^x w { + Zz wu l + + M^(^yWW{ + IjWl + <? 7 ^ + £yW>Vj) 

+ M^ 7 (CM^| + Zx^i + £ 7 W{ + + ^(lx + ^ + d)7{ ( A -2e) 

For linearization it is convenient to rewrite the last element as 
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(E^) 5 = [^(w 2 )j + $ 2 y{v\ + ^(H' 2 )j] +(jl + QltxtyMf + f x SJuw) s + ^»{] 

+ Y C&v 2 + + i]{u 2 + w 2 ), + £{1? + v 2 ),] + k{£ + + g)T t (A.2f) 

The elements of F Vl and have exactly the same form as those of E n , but with £ replaced by rj and C, 
respectively. 

The five elements of Ey 2 (again excluding the 1 jJRe r coefficient) are 


(E Ka >! — 0 (A.3a) 


(^V^2 ~ b b b b Cy v £ *b 

+ ^y(>ly^ v + *lx V r, + + ^JC V ^) + + VS, + ^Z U < + ^JC W ^ (A. 3b) 


(E[/ 2 )3 = 2 niy(rjyV v + CyV^) + + Vn + W, + ^ K C + + 

+ + f-* v n + + £jc v s) + MWn* v , + + C r v 5 + ^W f ) (A. 3c) 


(E^4 = z w v + tz w d + + vs? + ^ w c + £y v c + k w f) 

+ nixing + + ^ w c + £x w q) + a^(v,? + VS + £z v c + ^ w ;) ( AJd ) 


(E^) s = 2n[Z x (r] x UU^ + z x uu£ + iy(rtyW n + C^) + + C r w»v { )] 

+ + C y w c + C r iw c ) 

+ AtyiVxWn + + 'Jz VW »7 + + ^ W C + C Z VW Q ) 

+ ^(V% + V% + V™.? + + ^ w { + ^ w {) 

+ ntx (V - ^ + *»*% + r <z WU V + *fjc WM 'i| + + ^ w c + + Zx 1 ™# 

+ n^ y {fiyiM n + t] x uv n + r,^ + *\ y ww n + ZyliUt- + C x uv c + ZzWV/. + Z y ww£ 

+ ulA r, z uu n + Y} x uw n + ri z w n + r,yVW v + + t x UW Q -»- + C y vw^) 

+ k{Z x n x + fa + ^)T„ + KUx + Z& + (A-3e) 

A A 

The elements of FV 2 have exactly the same form as those of Ek 2 , but with £ replaced by r\ and rj replaced 

by <£. Similarly, the elements of Gk 2 have exactly the same form as those of Ek 2 , but with £ replaced by C 
and £ replaced by £. 
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